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PREFACE 

DTNSRDC  is  developing  a programming  language  for  describing 
structural  geometry.  This  language,  called  G-PRIME,  is  designed  to  assist 
scientists  and  engineers  in  translating  from  concepts  and  blueprints  to 
the  data  format  required  by  mathematical  simulation  programs.  G-PRIME 
must  be  able  to  accurately  represent  many  geometric  forms,  and  to 
reproduce,  scale,  move,  truncate,  combine,  and  find  intersections  of  these 
various  forms.  One  feature  which  greatly  simplifies  the  design  of 
G-PRIME,  and  one  that  probably  makes  it  unique,  is  that  all  curve  and 
surface  geometry  is  represented  using  just  one  mathematical  form-- 
B- spline  functions. 

A number  of  projects  within  the  Navy  ?buld  benefit  from  a unified 
set  of  curve  and  surface  manipulation  subroutines,  particularly  projects 
in  the  hull  representation  and  data  fitting  areas.  This  possibility 
prompted  us  to  isolate  G-PRIME' s basic  B-spline  manipulation  routines  and 
to  make  them  available  as  a separate  subroutine  library.  This  report 
contains  a discussion  of  the  mathematics  involved  in  the  use  of  B-spline 
functions,  a description  of  the  subroutines  contained  in  this  library,  and 
a sample  application  which  illustrates  a typical  use  of  the  subroutines. 

We  have  been  pleased  with  the  performance  of  these  subroutines  in 
our  applications.  It  should  be  noted,  however,  that  only  a first 
increment  of  the  intersection  capability  has  been  included  at  this  time. 
Further,  we  now  see  areas  in  which  improvements  can  be  made  and  in  which 
additional  capability  is  necessary.  This  report  concludes  with  a summary 
of  our  current  work  in  the  development  of  this  basic  B-spline  capability 
and  our  recommendations  for  future  directions. 
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ABSTRACT 

This  report  describes  a library  of  mathematical 
subroutines  for  defining  and  operating  on  B-spline 
curves  and  surfaces.  This  library  contains  sub- 
routines for  evaluating  B-spline  functions,  for 
using  B-splines  as  a basis  for  fitting  curves  and 
surface  data,  and  for  finding  the  intersections  of 
B-spline  curves  and  surfaces.  An  explanation  of 
the  general  theory,  a summary  of  important  properties, 
and  detailed  operating  instructions  for  each  sub- 
routine are  given.  A program  illustrating  a typical 
use  of  the  basic  subroutines  has  been  included  along 
with  computer -generated  plots  of  B-spline  surfaces 
and  intersection  curves. 

Although  the  basic  evaluation  subroutines  use 
previously  published  techniques,  the  fitting  and 
intersection  procedures  represent  effective  new 
approaches  to  the  treatment  of  these  old  problems. 


PARAMETRIC  B-SPLINE  CURVES  AND  SURFACES 

their  introduction  by  I.J.  Schoenberg  in  1946/  mathematical 
spline  techniques  have  become  increasingly  popular  for  smoothing  and 
interpolating  data  and  for  functional  approximation.  Much  of  the  popu- 
larity of  these  techniques  can  be  attributed  to  the  fact  that  mathemat- 
ically smooth  and  aesthetically  pleasing  shapes  are  easily  produced  using 
splines.  R.  Riesenfeld  has  proposed  using  spline  functions  for  curve  and 
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Schoenberg,  I.J.,  "Contributions  to  the  Problem  of  Approximations  of 
Equidistant  Data  by  Analytic  Functions,"  Quart.  Appl.  Math.,  vol.  4 (1946), 
pp.  45-99;  112-141. 
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surface  approximation  in  computer-aided  design  in  his  1973  thesis. 

The  design  method  which  he  proposes  is  patterned  after  Bezier's 
technique,'^  but  uses  a well-behaved  set  of  spline  functions  as  the  basis 
for  approximation  rather  than  the  Bernstein  polynomial  basis  usually 
associated  with  Bezier's  method.  The  advantages  of  using  these  basic 
splines  (or  B-splines)  include  numerical  stability,  computational 
efficiency,  and  the  ability  to  represent  a variety  of  common  curves  and 
surfaces  which  are  difficult  to  approximate  using  earlier  techniques. 
Riesenfeld's  thesis  provides  a good  introduction  to  both  Bezier's  method 
and  the  theory  of  B-splines. 

Bezier  and  Riesenfeld  choose  to  represent  curves  and  surfaces  by 
vector -valued  parametric  functions  instead  of  the  scalar  equations  which 
are  most  often  used  for  classical  curves  and  surfaces.  Although  it  is 
not  necessary  to  use  parametric  equations  with  B-splines,  the  parametric 
form  does  enlarge  the  family  of  admissible  shapes  and  makes  possible 
several  computational  short  cuts. 

Our  objective  has  been  to  produce  a library  of  FORTRAN  subroutines 
which  can  be  used  to  perform  most  elementary  operations  involving  B-spline 

curves  and  surfaces.  The  subroutines  described  in  this  report  include 

? 

(1)  an  implementation  of  the  techniques  described  by  Riesenfeld""  for 
evaluating  cubic  B-splines  having  a uniform  knot  vector,  (2)  routines  for 
fitting  data  with  B-spline  curves  and  surfaces,  and  (3)  routines  for 
finding  the  intersections  of  B-spline  curves  and  surfaces. 

For  the  reader  who  is  not  familiar  with  B-splines  the  following 
examples  may  help  to  illustrate  the  utility  of  such  a library. 

• A designer  wishes  to  produce  a transparent  plastic  canopy 
to  enclose  an  instrument  package.  Since  B-spline  functions  can  be  used 
directly  to  generate  a smooth  surface,  the  designer  needs  only  to  supply 
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Riesenfeld,  R. , "Application  of  B-Spline  Approximation  to  Geometric 
Problems  of  Computer-Aided  Design,"  University  of  Utah  Computer  Science 
Report  UTF.C-CSC-73-126. 

^ Bezier,  P. , Numerical  Control  - Mathematics  and  Applications  (translated 
by  A.R.  Forrest)-!  London:  .John  Wiley  and-  Sons,  1972. 
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a rectangular  mesh  of  points  which  describes  the  essential  shape  of  the 
canopy  and  its  boundaries  in  order  to  produce  a trial  surface.  The 
evaluation  routines  can  be  invoked  to  calculate  a number  of  points  on  the 
B-spline  surface  induced  by  the  data  given  in  the  mesh.  This  surface  may 
be  graphically  displayed  for  visual  evaluation,  it  may  be  used  with 
numerical  integration  techniques  to  calculate  the  volume  enclosed,  or  it 
may  be  used  to  check  clearances  for  the  instrument  package.  If  the  design 
is  not  satisfactory,  the  designer  can  change  it  by  modifying  the  points 
in  the  mesh.  When  a satisfactory  design  has  been  found,  the  evaluation 
routines  can  be  called  upon  to  supply  enough  points  on  the  surface  to 
produce  a mold  for  manufacturing  the  canopy. 

• An  engineer  is  laying  out  detailed  drawings  for  the 
fabrication  of  a plenum  which  involves  circular  pipes  and  an  oval  chamber 
meeting  at  oblique  angles.  He  needs  to  know  the  curves  of  intersection 
of  these  parts  for  cutting  the  material.  Points  on  the  curves  of  inter- 
section can  be  calculated  using  B-spline  functions  in  a two-step  process. 
First,  the  engineer  provides  data  points  on  each  of  the  intersecting 
pieces.  These  points  are  used  by  the  fitting  routines  to  determine 
B-spline  surfaces  which  represent  each  of  the  original  surfaces.  The 
intersection  subroutines  will  determine  points  on  the  curves  of  inter- 
section of  these  B-spline  surfaces. 

• A technician  would  like  to  have  experimental  data  plotted 
and  have  smooth  curves  drawn  which  interpolate  the  data  points  or  pass 
close  to  the  data  points.  The  fitting  routines  can  be  used  to  determine  a 
B-spline  curve  which  either  interpolates  the  data  or  best  fits  the  data  in 
the  least -squares  sense.  The  evaluation  routines  can  then  be  called  to 
obtain  sufficient  points  on  the  curve  for  smooth  plotting. 

The  remainder  of  this  section  describes  the  various  mathematical 
procedures  currently  available  in  this  library'. 
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EVALUATION  PROCEDURES 
The  B- Spline  Curve 

A parametric  B- spline  curve,  Q^P'.s] , is  said  to  be  induced  by  the 
polygon  P = P^.P^,...,?  m P°^n^s»  where  each  point  is  a vector  with  as 

many  components  as  the  dimensionality  of  the  space.  For  a given  polygon 
P,  and  parameter  value  s,  the  point  on  the  induced  B-spline  curve  is 
given  by 

m 

(yp;s]  = Z NM(s)P.  0 < s < 1 (1) 


where  N.  . (s)  is  the  i*h  normalized  cubic  B-spline  basis  function. 

1,4  4 

N.  .(si  is  computed  using  de  Boor's  recursive  procedure  : 

^ 9 ^ 


For  M-l 


(S) 


1 for  i 5 s- (m-l)<  i+1 
0 otherwise 


then  for  M = 2,3,4 


(2) 


Nf ,mCs)  = [(z-i)-NijM1(s)+(m+i-z)-N 


i+1, M-l 


(s) ] / (m-l) 


where  z = s-(m-l) 


For  closed  curves  a periodic  basis  is  required  which  can  be  obtained 
by  modifying  the  above  procedure  so  that  all  differences  and  indices  are 
computed  modulo  (m-l).  Whenever  0 ^ s < 1/ (m-l)  or  (m-l) -l/(m-l)  < s ^ 1 
the  basic  functions  may  be  modified  such  that  the  end  points  of  the 
polygon  P are  interpolated  by  the  B-spline  curve  and  the  slopes  at  the 
ends  of  the  curve  are  the  same  as  the  slopes  of  the  beginning  and  ending 
legs  of  the  polygon.  This  modification  is  equivalent  to  introducing  a 
non-uniform  knot  vector  in  the  spline  definition.  For  closed  curves, 
values  of  s outside  the  interval  0 < s - 1 are  replaced  by  s modulo  1.  For 
convenience,  values  of  s outside  the  interval  0 5 s 5 1 are  also  permitted 
for  open  curves.  In  this  event  basis  functions  are  constructed  to 
produce  a straight  line  extension  of  the  corresponding  leg  of  the  inducing 
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de  Boor,  C. , "On  Calculating  with  B-Splines,"  J.  Approx.  Theorv,  vol . 6 
( 1972 ),  pp.  50-62. 
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polygon.  (^[P-,0]  and  [ P ; 1 ] are  still  referred  to  as  the  "ends"  of  the 

curve,  however. 


A Linearly  Interpolated  Approximation  to  a B-Spline  Curve 

Frequently  it  is  necessary  to  evaluate  a B-spline  curve  for  some  set 

of  uniformly  spaced  parameter  values,  say  (s^.s^ s^) , and  Equation  (1) 

can  be  written  in  matrix  form: 


Q = [S]P  (3) 

where 

SiJ  • P-  [p1-P2 Pm'T-  ^ 

9 = [Qm[P:s1].  Qm[P;s2] ^[P;sk]1 


If  a matrix  of  B-snline  coefficients,  [S] , has  been  previously 

computed  and  [ P ; ] is  required  where  s2  - s|  - s^,  Q,  an  approximation 

to  Q may  suffice  for  certain  applications.  We  can  obtain  one  such 

approximation  by  linear  interpolation  of  the  coefficients  available  in 

[S] . If  the  two  parameter  values  used  in  computing  [SI,  which  bracket  s^, 

are  s . and  s . , then 

J 1+1 


Substituting  Equation  (1)  in  the  last  equation  gives 


or 

= £ NC)4(Si)Pc  (4 

The  B-Spline  Surface 

A parametric  B-spline  surface  is  said  to  be  induced  by  an  m by  n 
rectangular  mesh  of  points  [P] , where  each  P^j  is  a vector  with  as  many 
components  as  the  dimensionality  of  the  space. 

A point  on  a cubic  B-spline  surface  corresponding  to  the  pair  of 
parameter  values  (s,t)  is  given  by 
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(5) 


VIIP,;s,tl 


m n 

l E N (s)N.  . (t)P. . ; 0<S<1,  0<t<l 
i=l  j=l  1,4 


This  surface  is  a Cartesian  product  generalization  of  a B-spline  curve  and 
all  qualifications  imposed  on  the  curve  definition  apply  to  the  surface 
definition  as  well.  The  curves  Q^t  [P]  ;0 ,tl , [P]  ;1  ,t] , Q^t  [P]  ;s  ,0] , 

and  Q^[[P];s,l]  are  referred  to  as  the  "edges"  of  the  surface. 

To  evaluate  a B-spline  surface  for  a set  of  points  corresponding  to 
the  various  combinations  of  the  parameters  {s^ ^ , • . . ,s^}  and  *^2’ ‘ 
Equation  (5)  can  be  written  in  matrix  form: 

[0]  = [S]  [P]  [TT]  (6) 


where  S. . = N.  .(s.),  T..  = N.  .(t.),  [P]  is  the  inducing  mesh  of  points, 
1J  J 1 1J  J 1 

and  [C1  is  the  mesh  of  points  on  the  B-spline  surface. 

If  the  B-spline  coefficient  matrices  [S]  and  [T]  are  available,  a 

linearly  interpolated  approximation  to  the  B-spline  surface  can  be 

obtained  by  substituting  N.  , (z)  of  Equation  (4)  for  N-  . (z)  in  Equations 

J d J 

(5)  and  (6) . 

FITTING  PROCEDURES 


The  B-spline  approximation  to  a curve  or  surface  is  particularly 
useful  in  the  process  of  designing  objects  for  which  there  are  few 
quantitative  constraints  on  the  final  form  and  for  which  the  main  concerns 
are  gross  dimensions,  smoothness,  and  aesthetic  acceptability.  Such  an 
approximation  is  usually  not  satisfactory  as  a mathematical  representation 
of  the  shape  of  an  existing  object.  Many  of  the  advantages  of  B-spline 
techniques  can  still  be  brought  to  bear  on  this  problem  by  employing  least - 
squares  fitting  procedures  using  B-spline  functions  as  the  basis  for 
fitting  the  given  shape. 


Curve  Fitting 

The  problem  of  fitting  an  ordered  set  of  data  points,  Q,  with  a 
B-spline  curve  is  one  of  determining  the  inducing  polygon,  P,  in 

Q = [SIP 

such  that  the  error, 

e = |Q-Q| 
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2 

is  minimized  where  the  norm  |x|  = z x.  , and  [S]  is  the  B-spline 

all  i 1 

coefficient  matrix  used  in  Equation  (3). 

If  we  stipulate  that  there  are  at  least  as  many  points  in  Q as 
there  are  in  P,  we  can  employ  standard  linear  least -squares  fitting 
techniques  and  write  the  following  expression  for  P: 

P = [U]Q  (7) 

where 

[U]  = [ST  S]  'V  (8) 

Surface  Fitting 

The  generalization  of  a least -squares  fitted  curve  to  a surface  is 
more  easily  accomplished  by  introducing  operator  notation  and  zero 
indexing  for  the  sums.  If  the  cubic  B-spline  approximation  is  expressed 
as  a linear  operator,  (^,  on  C^[0,1],  where  d is  a positive  integer, 
and  the  curve  to  be  approximated  is  the  function  f(s),  the  approximation 
can  be  written  as 

m 

0 [f;s]  = l N.  .(s)  f (i/m) 

™ i=0  ’ 


The  least-squares  transformation  [U]  from  Equation  (8)  can  be  used 
to  define  a least-squares  fit  operator  F , such  that 

m k 

F [f;s]  = Z i N.  (s)U.  fU/k)  (9) 

m i=0  £=0  15 

If  the  data  point  Qf  from  Q is  substituted  for  each  f(E;/k),  Equation  (9) 
becomes 

m 

FmtQ;s]  = 2 4(s)Pi 
i=0  ’ 


P = [U]Q 

With  the  same  notation  the  Cartesian  product  bivariate  approximation 
of  a function  f(s,t)  is 

m n 

Vn[f;$’t]  = QmQn[f;s’tl  = ^ Ni,4(s)Nj  , 4 Ct ) f ( i/m , j /rO 
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The  Cartesian  product  generalization  of  the  least-squares  fitted  curve  is 
then 


F 

m,n 


[ f ; s , 1 1 

n k 


z z 

i=0  £=0 


n 8. 


I Z N 4(s)N  < 

j=0  n=0  1,4  J’ 


(tlUiCWjnf(f;/k’n/f) 


(10) 


If  the  data  point  Qr ^ from  [Q] 
Equation  (10)  becomes 


is  substituted  for  each  f U/k,nA) , 


m n 

= Z Z 

i=0  j=0 


Ni,4^Nj,4^Pij 


where 

fP]  « [U][Q][WT] 

[U]  = [sWi 

[W]  = tTTT] _1[TT] 


(11) 


(12) 


and  [S]  and  [T]  are  the  B-spline  coefficient  matrices  used  in  Equation  (6). 

For  certain  surface  design  applications  it  may  he  desirable  to  use 
fitting  techniques  for  one  parameter  and  a B-spline  approximation  for  the 
other.  This  can  be  accomplished  by  substituting  the  identity  matrix  for 
either  fU]  or  [IV]  in  Equation  (11),  depending  on  the  direction  in  which 
the  approximation  is  desired.  Once  [U]  and  [IV]  have  been  computed  for  a 
problem  of  a given  size,  the  fitting  procedure  for  each  set  of  data,  [Q]  , 
involves  only  matrix  multiplications.  Some  applications  can  exploit  this 
property  to  reduce  computing  time. 


INTERSECTION  PROCEDURES 


The  problem  of  finding  points  of  intersection  of  arbitrary  curves 
and  surfaces  is  essentially  that  of  solving  a system  of  nonlinear  equations. 
We  have  approached  the  task  of  finding  points  of  intersection  of 
parametric  cubic  R-spline  curves  and  surfaces  by  treating  them  as  arbitrary 
nonlinear  functions,  rather  than  taking  advantage  of  the  special 
properties  of  B-splines.  We  have  chosen  this  approach  for  two  reasons: 

(1)  the  B-spline  functions  used  for  geometric  applications  are  usually  very 
well  behaved  and  efficient  numerical  techniques  for  solving  systems  of 
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well-behaved  nonlinear  equations  are  readily  available;  and  (2)  one 
solution  procedure  can  be  used  for  finding  intersections  involving  curves, 
surfaces,  curves  with  surfaces,  and  for  finding  the  point  on  a curve  or 
surface  which  is  closest  to  a given  spatial  point.  We  feel  that  the 
straightforward  programing  that  can  be  achieved  using  this  approach  more 
than  offsets  any  computing  overhead  incurred  by  not  taking  advantage  of 
special  properties. 

After  evaluating  several  iterative,  nonlinear  solvers  we  have 
chosen  Brown's  algorithm^  for  its  efficiency  and  reliability.  This 
algorithm  employs  a nonlinear,  least-squares  minimization  technique  that 
does  not  require  explicit  derivatives  of  the  functions  to  be  minimized. 

The  Intersection  of  Two  Plane  Curves 

Although  a procedure  for  finding  the  points  of  intersection  of  two 
planar  curves  is  not  available  in  the  current  version  of  B-spline 
subroutine  library,  this  problem  illustrates  all  the  essential  features 
of  the  existing  intersection  procedures. 

Consider  the  two  parametric  curves 


w - 


xl(Si)| 

>,i(si)l 


and  ?2(s2) 


IX2(S2) 

ly2(S2) 


(13) 


shown  in  Figure  1.  The  point  of  intersection  may  be  found  by  finding  a 
pair  of  parameter  values,  (s^,s2),  that  makes  the  residual  function 


R = 


|x2(s2) 

|y2(s2) 


X^Sj) 

Yl(Si) 


f2(s2^  ' 


(14) 


zero.  R is  the  function  that  is  to  be  minimized  by  the  iterative  solver. 
For  R-spline  curves  the  f^(s)'s  are  calculated  using  one  of  the  B-spline 
evaluation  procedures. 

Iterative  procedures  require  some  initial  value  for  each  of  the 
independent  variables.  When  there  is  only  one  point  of  intersection,  as 


Brown,  K.M. , "Derivative  Free  Analogues  of  the  Levenberg-Marquardt  and 
Gauss  Algorithms  for  Nonlinear  Feast  Squares  Approximations,"  Numerische 
Mathematik , vol.  18,  pp.  289-297,  1972. 
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in  Figure  1,  almost  any  pair  of  initial  values,  which  are  noc  both  zero, 
will  yield  convergence  to  the  correct  solution. 

If  the  two  curves  have  no  point  of  intersection  within  the  range 
defined  by  parameter  values  between  zero  and  one  and  the  linear 
extensions  of  the  B-spline  curves  do  intersect  (see  Evaluation  Procedure 
section),  the  algorithm  will  quickly  converge  to  a point  outside  the 
range.  If  the  extended  curves  do  not  intersect,  the  minimization  algorithm 
will  choose  the  two  closest  points  on  the  curves  as  the  solution  and  will 
terminate  with  a non-zero  residual  function. 

Locating  the  Closest  Point  on  a Curve 

If  there  are  multiple  intersection  points,  as  in  Figure  2,  the  user 
must  select  between  them  by  specifying  initial  parameter  values  which  are 
closest  to  the  desired  intersection  point.  If  there  are  only  two  inter- 
section points,  this  parameter  estimate  may  not  be  an  unreasonable  request. 
Clearly,  a spatial  estimate  of  the  location  of  the  desired  intersection 
would  be  more  natural  for  the  user.  Further,  spatial  estimates  seem  to  be 
mandatory  if  this  capability  is  to  be  used  effectively  with  computer 
graphics  devices  that  have  light  pens  or  other  graphical  means  of  input. 

If  the  user  gives  a spatial  point,  x^  = (x0,Vq),  the  point  on  the 
curve,  fj(s),  that  is  closest  to  the  point  x^  can  be  found  by  determining 
the  parameter  value,  s^ , that  minimizes  the  residual  function 
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Figure  2 - Two  Curves  with  Two  Intersection  Points 


\ (sx) 

n(si) 


fj(Si)  - XQ 


Again  the  iterative  nonlinear  minimization  procedure  is  used  to 
solve  for  s^.  If  the  procedure  is  repeated  with  the  second  curve,  in  order 
to  determine  a parameter  estimate,  s2,  the  required  initial  parameter 
values  will  have  been  obtained  to  begin  to  solve  for  the  desired  inter- 
section point. 

The  Intersection  of  Three  Surfaces 

The  procedure  for  finding  the  point  of  intersection  of  three  surfaces 
is  the  same  as  that  described  for  two  plane  curves  except  that  a pair  of 
parameters,  (s^,t^),  must  be  determined  for  each  of  the  surface  functions, 
f i Csi » t i) , which  makes  residual  function 

x2^s2,t2^  * 

y2(S2’t2)  ' f (s  ,t2)-f  (s  t,) 

WV  ' zi(srV  2 2 2 111 

x3(s3’V  ‘ x2Cs2 »t2)  f rc  t w r ♦ i 

y3(s3,t3)  - y2(s2,t2)  i3lS3’V'i2lS2’V 

Z3(^S3,t3-)  ' z2^S2,t2-) 


( 


To  locate  the  point  on  a surface  f^(Sj.t^)  that  is  closest  to  a given 
spatial  point,  x^  = (x^.y^.z^),  a pair  of  parameters,  (s^.t^J,  must  be 
found  which  will  minimize  the  residual  function 


R = 


x,(s1  .tj) 

new 

zi(W 


■ fWh’^o 


(17) 


Limitations  of  the  Method  for  the  Three-Surface  Problem 

It  should  be  noted  that  with  this  method  of  finding  points  of  inter- 
section it  is  not  too  difficult  to  pose  a problem  that  will  fail.  The 
most  obvious  example  can  be  seen  in  Figure  3,  which  contains  two  plane 
curves  which  intersect  at  one  point  and  come  very’  close  to  intersecting  at 
another  (stationary  point).  If  the  initial  parameter  values  are  associated 
with  points  that  are  closer  to  the  stationary  point  than  the  actual  inter- 
section, the  method  will  usually  converge  to  the  stationary’  point  and 
terminate,  indicating  failure  to  find  an  intersection.  Problems  of  this 
type  can  be  avoided  by  treating  them  as  multiple  intersection  point 
problems.  Other  situations  are  possible  which  can  be  avoided  only  by 
subdividing  the  problem  to  restrict  the  region  of  definition. 


Figure  3 - An  Intersection  Point  and  a Stationary  Point 
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The  Intersection  of  Two  Surfaces 

The  intersection  capability  of  the  B-spline  subroutine  library 
contains  procedures  for  finding  points  on  the  curve  of  intersection  of 
two  B-spline  surfaces.  These  points  may  be  used  directly  or  a fitted 
B-spline  curve  may  be  obtained  using  the  procedures  described  earlier. 

The  two -surface  intersection  procedure  has  two  major  processing 
steps.  First,  an  edge  analysis  is  performed  to  determine  the  general 
location  and  extent  of  the  intersection  curve.  The  results  of  this 
analysis  are  in  terms  of  ranges  of  parameter  values  on  each  of  the 
surfaces.  The  second  step  involves  traversing  the  curve  by  incrementing 
a selected  parameter  over  the  range  of  the  curve  and  determining  a point 
of  intersection  for  each  value  of  that  parameter. 

Edge  Analysis 

To  find  a curve  of  intersection  of  the  two  surfaces  f^s^.tp  and 
f 2 (s2 ’ t2^  some  insight  into  the  nature  of  such  a curve  can  usually  be 
gained  by  determining  the  points  at  which  the  boundary  of  one  surface 
pierces  the  other.  To  accomplish  this  a procedure  similar  to  the  one 
described  for  finding  points  of  intersection  can  be  used.  Determining 
the  existence  of  these  piercing  points  requires  the  investigation  of  eight 
possible  combinations,  each  of  which  involves  an  edge  of  one  surface  and 
the  other  surface  in  its  entirety.  The  respective  residual  functions  for 
these  eight  combinations  are  as  follows: 

1*1  = f1(0,t1)-f2(s2,t2) 

?2  = -l^1’tl^*-2^s2,t2^ 

-3  = ?l(sl’°)'-2('S2’t2') 

R4  = f1(s1,l)-f2Cs2,t2) 

After  the  edge  analysis,  problems  with  two  distinct  edge  pierce  points 
(and  certain  closed  surface  problems  with  one  edge  pierce  point)  will 
be  admitted  for  further  processing.  Any  other  outcome  will  cause  the 
procedure  to  terminate,  indicating  failure. 


r5  = f1(s1,t1)-f2(0,t2) 

-6  = ~l^Sl,tP'-2(-1,t2^ 
1*7  = f1(s1,t1)-f2(s2,0) 

-8  = -l(sl,tl)‘-2(s2’1) 


(18) 
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Computing  Points  on  the  Curve  of  Intersection 

The  parameter  values  associated  with  the  two  edge  pierce  points  are 
examined  and  the  one  with  the  greatest  range  is  selected  to  control  the 
traversing  of  the  intersection  curve.  If  the  selected  parameter,  s^,  has 
the  range  a 1 S2  1 b and  the  user  has  requested  that  k+1  points  be  found 
along  the  curve,  the  points  are  determined  by  submitting  the  following 
sequence  of  residual  functions  to  the  iterative  solver: 

R = f x (sx ,t j)  - f2(c,t2)  (19) 

where  c = a+(b-a) -i/k,  i = 0,1,..., k 

When  one  or  three  edge  pierce  points  is  found  and  one  of  the  surfaces  is 
a closed  surface,  the  closed  surface  parameter  associated  with  the 
periodic  B-spline  functions  is  chosen  to  control  the  traversing  of  the  curve. 

Limitations  of  the  Method  for  the  Two-Surface  Problem 

Current  two-surface  intersection  procedures  admit  only  problems  with 
one  or  two  edge  pierce  points.  This  limitation  excludes  any  problem  with 
multiple  curves  of  intersection  and  those  problems  in  which  a region  of  one 
surface  bulges  through  the  other  surface.  These  problems  can  usually  be 
treated  by  restricting  the  regions  of  definition  or  by  subdividing  the 
problem. 

PROPERTIES  OF  B- SPLINE  FUNCTIONS 

The  application  of  B-spline  techniques  is  quite  straightforward  and 

extensive  mathematical  insight  is  not  required  for  their  use.  Of  course, 

an  awareness  of  the  major  properties  of  B-spline  functions  should  help  the 

user  avoid  difficulties  and  perhaps  open  new  avenues  for  their  use.  A 

short  summary  of  properties  is  given  below.  Most  of  the  proofs  have  been 

? 

sketched  by  Risenfeld“  or  can  be  obtained  by  inspection  of  the  B-spline 
definitions. 

(1)  If  the  evaluation  procedures  are  applied  to  data  points  taken 
from  a given  curve,  the  induced  B-spline  curve  will  approximate  that  given 
curve  with  a certain  degree  of  accuracy.  As  defined,  B-spline  functions 
yield  a convergent  approximation;  hence,  the  more  data  points  used,  the 
better  the  approximation.  The  cubic  B-spline  approximation  is  also  a 
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convergent  approximation  for  the  first  three  derivatives  of  a function. 

(2)  In  regions  of  high  curvature,  more  data  points  will  be  required 
to  obtain  a given  level  of  approximation  than  in  gently  curving  regions. 

(3)  A B-spline  curve  will  either  coincide  with  the  inducing  polygon 
or  lie  on  the  concave  side  of  that  polygon  (convex  hull  property). 

(4)  The  cubic  B-spline  curves  described  here  have  continuous  second 
derivatives  everywhere  and  continuous  third  derivatives  everywhere  with 

the  possible  exception  of  the  points  0,  l/(m-l),  2/(m-l) 1.  A 

repeating  data  point  in  the  polygon  can  reduce  the  differentiability  by 
one  at  that  point.  Cusps  or  sharp  comers  can  be  induced  in  a B-spline 
curve  by  specifying  the  same  data  point  three  times. 

(5)  If  three  points  of  a polygon  lie  in  a straight  line,  the  induced 
curve  will  pass  through  the  middle  point  exactly. 

(6)  It  is  apparent  that  can  be  considered  a linear  operator  since 


qjaP  + bR;s]  = a(^[P;s]  ♦ b^JR-.s] 

This  implies  that  one  can  apply  a linear  transformation  to  a polygon  and 
use  that  transformed  polygon  to  induce  a B-spline  curve  which  is  equivalent 
to  applying  the  linear  transformation  to  the  B-spline  curve  induced  by  the 
original  polygon.  Relevant  linear  transformations  include  rotation, 
translation,  scaling,  and  projection. 


Each  of  the  above  properties  can  be  extended  to  B-spline  surfaces 
using  the  Cartesian  product  generalization  of  the  curve  definition. 


! 


BASIC  B-SPLINE  SUBROUTINES 


Our  implementation  of  the  B-spline  manipulation  procedures  is  in  the 
form  of  a library  of  FORTRAN  subroutines.  The  order  of  the  previous  section 
has  been  followed  to  group  the  subroutine  descriptions  according  to  the 
functions:  Evaluation,  Fitting,  and  Intersection.  A fourth  category, 
Utility,  describes  subordinate  routines  which  are  used  only  indirectly  in 
the  B-spline  manipulation  processes. 
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Although  there  are  three  primary  evaluation  subroutines,  the  main 
most  general  one  for  both  curves  and  surfaces  is  BSEVL1 . The  other 
routines,  BSEVL2  and  BSEVL3,  are  somewhat  more  efficient  under  special 
conditions.  Applications  which  involve  several  surfaces,  each  defined 
by  the  same  size  mesh,  and  which  are  to  be  evaluated  for  a fixed  set  of 
parameter  values  may  substitute  subroutine  BSLVL2  for  BSEVL1  to  reduce 
computation.  Subroutine  BSEVL3  could  also  be  substituted  for  BSIiVLl , 
but  since  BSEVL3  only  computes  approximations  to  B-spline  functions 
there  are  very  few  applications  which  can  benefit  from  the  use  of  this 
subroutine.  The  one  fitting  subroutine  included  in  the  library,  FITBS2, 
is  capable  of  performing  all  the  fitting  procedures  described  in  the 
previous  section  for  both  curves  and  surfaces. 

Only  two  intersection  subroutines  have  been  included  at  this  time. 

One,  INT3S,  finds  a point  of  intersection  of  three  B-spline  surfaces.  The 
other,  1NT2S,  finds  points  on  the  curve  of  intersection  of  two  B-spline 
surfaces.  INT3S  is  quite  general,  solves  most  problems  of  practical 
interest,  and  handles  many  pathological  situations  as  well.  INT2S  is 
restricted  to  problems  which  have  one  continuous  curve  of  intersection 
that  includes  at  least  one  point  on  one  of  the  edges  of  the  two  intersecting 
surfaces.  This  subroutine  will  handle  a good  range  of  practical  problems, 
and  a much  broader  class  of  problems  can  be  treated  if  the  user  is  willing 
to  subdivide  to  meet  the  stated  restrictions. 

Subroutines  INT2S  and  INT3S  can  also  be  accessed  with  simplified 
calling  sequences  by  using  subroutines  INT2SX  and  INT3SX,  respectively. 

The  utility  subroutines  may  be  used  independently  to  perform  their  various 
functions,  but  in  general  they  have  been  programmed  to  be  efficient  for  the 
tasks  at  hand  and  arc  not  necessarily  the  type  one  would  find  in  a general 
mathematical  subroutine  library. 

GLOBAL  CONVENTIONS 

All  the  B-spline  mathematical  subroutines  described  in  this  report 
adhere  to  the  following  conventions: 

fl)  Each  rectangular  mesh  of  points,  that  contains  data  points  taken  from 
a surface,  or  points  used  to  induce  a B-spline  surface,  must  be  stored  in 
a FORTRAN  array  dimensioned  as  follows: 
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P(MD,  Ml,  Nl) 

where  MD  is  the  dimensionality  of  the  space,  Ml  is  the  number  of  mesh 
points  along  an  edge  of  the  mesh  associated  with  the  variation  of  the 
first  parameter,  and  Nl  is  the  number  of  mesh  points  along  an  edge  of  the 
mesh  associated  with  the  variation  of  the  second  parameter.  The  first 
parameter  is  referred  to  as  the  "s"  parameter  and  the  second  parameter  is 
referred  to  as  the  "t"  parameter  (see  Figure  4) . 

\ 
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Figure  4 - Surface  Mesh  Conventions 


(2)  Each  polygon  of  points,  that  contains  data  points  taken  from  a curve, 
or  points  used  to  induce  a B-spline  curve,  must  be  stored  in  a FORTRAN 
array  dimensioned  as  follows: 

P(MD,  Ml) 

where  MD  is  the  dimensionality  of  the  space,  and  Ml  is  the  number  of  points 
in  the  polygon.  Note  that  this  polygon  is  equivalent  to  a surface  mesh 
with  either  Ml  = 1 or  Nl  = 1. 

(3)  An  integer  code  word,  ITYPF,  is  associated  with  each  mesh  or  polygon 
of  points  used  to  define  a B-spline  curve  or  surface.  The  type  of  B-spline 


curve  or  surface  desired  is  indicated  by  setting  ITYPE  according  to  the 
values  in  Table  1. 


TABLE  1 - CURVE  AND  SURFACE  TYPE.  CODES 


I TYPE 

CURVE.  OR  SURFACE  TYPE 

1 

Open  Curve 

2 

Closed  Curve 

3 

Open  Surface 

4 

s -Closed  Surface 

5 

t -Closed  Surface 

6 

s-  and  t -Closed  Surface 

(4)  A polygon  used  to  define  a closed  curve  must  have  both  points  of 
closure  given,  even  though  the  second  point  is  redundant.  A mesh  used  to 
define  a closed  surface  must  also  have  the  points  of  closure  repeated. 

For  a given  direction  a surface  will  be  either  open  or  closed.  No  hybrid 
combinations  can  be  represented. 

(5)  Any  real  number  is  valid  as  a parameter  value;  however,  for  non- 
periodic functions  a parameter  value  which  is  less  than  zero  or  greater 
than  one  yields  a point  on  a linear  extrapolation  of  the  function. 

(6)  Cubic  B-spline  coefficients  for  a given  parameter  value  are  normally 
packed  into  a four -word  array,  CN,  and  a packing  index,  LET,  is  associated 
with  those  coefficients.  Details  of  the  packing  algorithm  are  included 
with  the  description  of  subroutine  MKSPLN. 

For  most  applications  only  a few  of  the  library  subroutines  will  be 
called  directly  by  the  user.  For  example,  in  a program  involving  surfaces 
the  user  is  required  to  read  or  generate  the  mesh  of  points,  [Pj  , to 
define  each  surface.  He  may  then  determine  a new  mesh  of  points  which  will 
induce  a B-spline  surface  that  best  fits  his  original  set  of  points  by 
calling  subroutine  FITBS2.  Surface  meshes,  regardless  of  their  source,  arc 
all  used  in  the  same  way  in  the  application  portion  of  the  program.  There, 
subroutine  BSEVL1  can  be  called  to  obtain  the  coordinates  of  points  on  the 
surface  and  the  intersection  subroutines  can  be  used  to  find  intersection 
points  and  curves. 
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EVALUATION  SUBROUTINES 
Subroutine  Name 


Page 


BSCMAT 

20 

BSEVL1* 

21 

BSEVL2* 

22 

BSEVL3* 

24 

MKSPLN 

26 

PLSPLN 

29 

BSCMAT  - Subroutine  Description 


Function:  To  compute  a cubic  B-spline  coefficient  matrix. 

lintry  Point:  BSCMAT 


Calling  Sequence:  CALL  BSCMAT (CS ,M1 ,K1 , ICLOSL) 

Calling  Arguments: 

Name  (Attributes! 

CS  (mixed/array/output) 


Ml  (integer/ input) 

K1  (integer/ input) 
ICLOSL  (integer/ input) 


Contents 

Cubic  B-spline  coefficient  matrix, 
[S] , in  packed  form.  Dimensioned 
CS(5,K1) . 

Row  dimension  of  full  [S]  matrix. 
Column  dimension  of  [S]  matrix. 
Flag  word.  A value  of  zero 
indicates  that  a non-periodic 
B-spline  basis  is  to  be  used  to 
compute  coefficients.  A value  of 
one  indicates  that  a periodic 
basis  is  to  be  used. 


COMMON  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Sub rout ines : MKSPLN 

Method:  Subroutine  MKSPLN  is  called  to  compute  rows  of  the  coefficient 

matrix  for  the  K1  parameter  values  0,  1 . / (kl - 1) » 2 . / ( K1 - 1) , . . . ,1 . The  non- 
zero elements  of  each  row  are  stored  in  the  first  four  words  of  each  row 
of  CS  and  the  associated  integer  packing  index  is  stored  in  the  fifth  word 
of  each  row. 

Remarks : 

(1)  See  description  of  subroutine  MKSPLN  for  matrix  packing 
algorithm. 
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BSEVL1  - Subroutine  Description 


Ev; 


Function:  To  evaluate  a parametric  cubic  B-spline  curve  or  surface 

function  for  a given  set  of  parameter  values,  yielding  the  coordinates  Of 
a point  on  the  function. 

Fjitry  Point:  BSEVL1 

Calling  Sequence:  CALL  BSEVL1 (P,MD, Ml ,N1 ,S,T, ITYPE, XYZ) 

Calling  Arguments: 


Name  (Attributes) 

P (real/array/ input) 

MD,M1,N1  (integer/ input) 
S (real /input) 


T (real /input) 


ITYPE  (integer/ input) 


Contents 

Polygon  or  mesh  defining  a B- 
spline  curve  or  surface. 
Dimensions  of  array  P. 

Parameter  value  of  point  to  be 
evaluated  on  B-spline  curve  or 
first  parameter  of  point  on  a 
surface. 

Second  parameter  value  of  point 
to  be  evaluated  on  B-spline 
surface. 

Code  word  indicating  type  of 
B-spline  curve  or  surface  (see 
Table  1). 

Coordinates  of  point  evaluated. 
Dimens ioned  XYZ (MD) . 


XYZ  (real /array/output) 

COf-MON  Areas : None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  MKSPLN 

Method:  This  subroutine  evaluates  Equation  (1)  if  a polygon  defining 

a B-spline  curve  has  been  given  or  Equation  (2)  if  a mesh  defining  a B- 

spline  surface  has  been  given.  Subroutine  MKSPLN  is  called  to  compute  the 

B-spline  coefficients,  N.  .(s)  and  N.  . (t). 

1 ^ J 

Remarks : 

(1)  The  dimension  N1  and  the  parameter  value  T are  not  referenced 
if  a curve  is  indicated  by  ITYPE. 
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BSEVL2  - Subroutine  Description 

Function:  To  evaluate  a parametric  cubic  B-spline  curve  or  surface 

function  using  precomputed  coefficient  matrices,  yielding  the  coordinates 
of  a point  on  the  function. 

Entry  Point : BSEVL2 

Calling  Sequence:  CALL  BSEVL2(P,MD,M1 ,N1 ,K,L,CS,LS1 ,CT,LT1 , 

ITYPE.XYZ) 


Name  (Attributes) 

P (real/array/ input) 

MD ,M1  ,N1  ( integer/ input ) 
K (integer/ input) 


L (integer/ input) 


CS  (mixed/array/ input) 


LSI  (integer/ input) 


CT  (mixed/array/input) 


LT1  (integer/ input) 


Contents 

Polygon  or  mesh  defining  a B- 
spline  curve  or  surface. 

Dimensions  of  array  P. 

Index  of  the  parameter  value  s^ 
for  which  the  function  is  to  be 
evaluated. 

Index  of  the  parameter  value  t^ 
for  which  the  function  is  to  be 
evaluated. 

B-spline  coefficient  matrix  [S] 
of  Equations  (3)  and  (6)  in  packed 
form.  Dimensioned  CS(5,LS1). 

Column  dimension  of  [S]  and  the 
number  of  uniformly  spaced 
parameter  values  s^  for  which  the 
B-spline  coefficients  have  been 
calculated . 

B-spline  coefficient  matrix  [T] 
of  Equation  (6)  in  packed  form. 
Dimensioned  CT(S,LT1). 

Column  dimension  of  [T]  and  the 
number  of  uniformly  spaced  parameter 
values,  tj,  for  which  the  B-spline 
coefficients  have  been  calculated. 
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ITYPE  (integer/ input)  Code  word  indicating  type  of 

B-spline  curve  or  surface  (see 
Table  1) . 

XYZ  (real/array/output)  Coordinates  of  point  evaluated. 

Dimensioned  XYZ(MD). 

COMMON  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  None 

Method:  This  subroutine  evaluates  Equation  (1)  if  a polygon  defining 

a B-spline  curve  has  been  given  or  Equation  (5)  if  a mesh  defining  a B- 
spline  surface  has  been  given.  The  coefficients  N^  ^is^)  and  N^  ^ ( t „ ) are 
obtained  from  the  arrays  CS  and  CT  respectively. 

Remarks : 

(1)  The  dimension  Nl,  the  parameter  index  L,  the  coefficient 
matrix  CT,  and  the  dimension  LT1  are  not  referenced  if  a curve  is  indicated 
by  ITYPE. 

(2)  The  coefficient  matrices  CS  and  CT  can  be  computed  by 
subroutine  BSCMAT. 


} 


V 
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BSEVL3  - Subroutine  Description 


Function:  To  evaluate  a linearly  interpolated  approximation  to  a 

cubic  B-spline  curve  or  surface  function,  yielding  the  coordinates  of  a 
point  on  the  approximate  function. 

Entry  Point:  BSEVL3 

Calling  Sequence:  CALL  BSFVL3(P,MD,M1 ,N1 ,S,T,CS,LS1 ,CT,LT1 , 

ITYPE,XYZ) 


Calling  Arguments: 

Name  (Attributes) 

P (real/array/ input) 

MD,M1,N1  (integer/ input) 
S (real/input) 


T (real/ input) 


CS  (mixed/array/ input) 


LSI  (integer/input) 


CT  (mixed/ array/ input) 


LT1  (integer /input) 


Contents 

Polygon  or  mesh  defining  a B-spline 
curve  or  surface. 

Dimensions  of  array  P. 

Parameter  value  of  point  to  be 
evaluated  on  the  approximate  curve 
or  first  parameter  of  point  on  a 
surface . 

Second  parameter  value  of  point  to 
be  evaluated  on  approximate 
surface. 

B-spline  coefficient  matrix  [S] 
of  Equations  (3)  and  (6)  in  packed 
form.  Dimensioned  CS(5,LS1). 

Column  dimension  of  [Si  and  the 
number  of  uniformly  spaced 
parameter  values  s ^ for  which  the 
B-spline  coefficients  have  been 
calculated. 

B-spline  coefficient  matrix  [T]  of 
Equation  (6)  in  packed  form. 
Dimensioned  CT(S,LT1). 

Column  dimension  of  [T]  and  the 
number  of  uniformly  spaced  parameter 
values,  t • , for  which  the  B-spline 
coefficients  have  been  calculated. 
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Table  1). 

XYZ  (real/array/output)  Coordinates  of  point  evaluated. 

Dimensioned  XYZ(MD). 

CCMMON  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  PLSPLN 

Method : This  subroutine  evaluates  Equation  (4)  if  a polygon  defining 

a B-spline  curve  has  been  given  or  a variant  of  Equation  (5)  if  a mesh 

defining  a B-spline  surface  has  been  given.  Subroutine  PLSPLN  is  called 

to  compute  the  approximate  B-spline  coefficients,  N.  . (s)  and  N.  . (t). 

i ***  J i’ 

Remarks : 

(1)  The  dimension  N1 , the  parameter  T,  the  coefficient  matrix  CT, 
and  the  dimension  LT1  will  not  be  referenced  if  a curve  is  indicated  by 

I TYPE. 

(2)  The  coefficient  matrices  CS  and  CT  can  be  computed  by  sub- 
routine BSCMAT. 

(3)  This  subroutine  was  designed  to  quickly  calculate  an 
approximation  to  a B-spline  function  for  the  early  stages  of  iterative  non- 
linear solution  routines.  In  tests,  it  has  executed  at  about  twice  the 
speed  of  subroutine  BSEVL1,  which  was  not  fast  enough  to  be  beneficial  to 
the  overall  process. 


MKSPLN  - Subroutine  Description 


To  compute  B-spline  coefficients  for  a given  parameter 


CALL  MKSPLN (S ,CN ,M1 ,MBORD , ICLOSE ,LFT) 


Function: 
value. 

Entry  Point:  MKSPLN 

Calling  Sequence: 

Calling  Arguments: 

Name  (Attributes) 

S (real/ input) 

CN  (real /array/output) 

Ml  (integer/input) 

MBORD  (integer/ input) 

ICLOSE  (integer/ input) 


LFT  (integer /out put) 


Contents 

Parameter  value  for  which 
coefficients  are  to  be  computed. 
B-spline  coefficients  in  packed 
form.  Dimensioned  CN (MBORD) . 

Length  of  full  (unpacked)  vector 
of  coefficients. 

Order  of  B-spline  basis  function 
to  be  used  to  compute  coefficients. 
Flag  word.  If  zero,  a non- 
periodic or  open  B-spline  basis  is 
to  be  used.  If  one,  a periodic 
or  closed  B-spline  basis  is  to  be 
used, 

A matrix  packing  index.  See 
Method  for  description. 


COMON  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  None 

Method:  Non-zero  B-spline  coefficients  are  computed  using  Formula  (2) 

and  stored  in  the  array  CN.  A matrix  packing  index,  LIT,  is  computed  for 
later  reconstruction  of  a full  vector.  LFT  = S*FLOAT (Ml - 1 ) + 1 . for  a 
non-periodic  basis.  For  a periodic  basis  the  computation  of  LFT  is  made 
modulo  (Ml-1). 

Remarks : 

(1)  For  cubic  B-spline  functions  MBORD  must  be  set  to  4. 

(2)  A specific  B-spline  coefficient  used  in  Formula  (2)  can  be 
obtained  from  the  array  CN  using  the  following  algorithm  (note  that  the 
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argument  MBORD  is  the  same  as  the  subscript  M in  the  formula) : 

a.  For  a non-periodic  basis  the  coefficient  M(S)  is  stored  in 
CN(IPOINT)  where  IPOINT  = i-LFT+(MB0RD+l)/2  whenever 

1 < IPOINT  < MBORD.  Otherwise,  Ni  (S)  = 0. 

b.  For  a periodic  basis  the  coefficient  ^(S)  is  stored  in 
CN(IPOINT)  where  IPOINT  is  computed  as  in  the  non-periodic 
case,  except  the  computation  is  performed  modulo  (Ml-1).  Note 
that  1 5 i . (Ml-1)  for  a periodic  basis. 

If  the  index  i is  required,  it  may  be  obtained  from  the  index  IPOINT  and 
the  parameter  value  using  the  formula: 

i = IPOINT  + LFT  - (MBORD+D/2 

For  a non-periodic  basis  values  of  i outside  the  interval  [1,M1]  must  be 
excluded.  For  a periodic  basis  i must  be  evaluated  modulo  (Ml-1)  when 
i > (Ml-1)  and  taken  as  i modulo  (Ml-1)  + (Ml-1)  when  i < 1. 

(3)  Modified  coefficients  which  cause  non-periodic  B-spline  functions 
to  pass  through  the  end  points  of  a defining  polygon  are  calculated  for 
cubic  B-splines  (MBORD=4).  An  open  curve,  for  example,  to  be  represented 
by  data  points  associated  with  Ml  uniformly  spaced  parameter  values, 
requires  extra  data  near  the  end  points  for  a non-periodic  B-spline  basis. 
Specifically  required  are  data  at  points  associated  with  the  parameter 

Values  S = 1. /FLOAT (3.* (Ml-1))  and 

S = 1 . - 1 . /FLOAT ( 3 . * (Ml  - 1 ) 

Roughly,  this  represents  data  at  one  third  the  distance  between  the  end 
point  and  its  adjacent  point.  These  values  are  obtained  by  linear  inter- 
polation. This  method  produces  a B-spline  curve  which  (1)  interpolates 
the  end  points  of  the  inducing  polygon,  and  (2)  has  the  same  slope  at  the 
endpoints  as  does  the  end  leg  of  the  polygon. 

(4)  Although  we  chose  to  use  normalized  B-spline  functions  (curves 
and  surfaces  were  defined  only  for  parameter  values  between  zero  and  one) , 
it  is  sometimes  useful  to  have  smooth  continuation  of  a curve  or  surface 
outside  the  range  of  definition.  Subroutine  MKSPLN  does  provide  for  linear 
extrapolation  of  open  B-spline  forms  when  parameter  values  outside  the 
range  zero  to  one  are  given.  For  parameter  values  which  are  less  than  zero 
the  coefficients  are  packed  in  array  CN  as  they  would  be  for  a parameter 
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value  of  zero.  Similarly,  for  parameter  values  greater  than  one,  the 
coefficients  are  packed  as  they  would  be  for  a parameter  value  of  one. 
For  periodic  B-spline  functions,  the  parameter  value  used,  SS,  will  be 
computed  from  a given  value  S by  the  formula. 

SS  = AM0D(AM0D(S,1.)+1. ,1.) 
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PLSPLN  - Subroutine  Description 

Function:  To  compute  approximate  B-spline  coefficients  for  a given 

parameter  value  by  linearly  interpolating  from  precomputed  coefficients. 
Entry  Point:  PLSPLN 

Calling  Sequence:  CALL  PLSPLN (S,CN, Ml , MBORD, CS,K1 , ICLOSE.LFT) 


Name  (Attributes)  Contents 

S (real/ input)  Parameter  value  for  which 

coefficients  are  to  be  computed. 

CN  (real/array/output)  B-spline  coefficients  in  packed 

form.  Dimensioned  CNCMBORD+l). 

Ml  (integer/input)  Length  of  full  (unpacked)  vector 

of  coefficients  and  row  dimension 
of  [SI  matrix. 

MBORD  (integer/ input)  Order  of  B-spline  basis  function 

to  be  used  to  compute 
coefficients . 

CS  (mixed/array/ input)  B-splme  coefficient  matrix,  [S]  , 

in  packed  form.  Dimensioned 
CS(5,K1) . 

K1  (integer/ input)  Column  dimension  of  [S]  matrix. 

I CLOSE  (integer/ input)  Flag  word.  If  zero,  a non- 

periodic or  open  B-spline  basis 
is  to  be  used.  If  one,  a periodic 
or  closed  B-spline  basis  is  to  be 
used. 

COMMON  Areas:  None 

FORTRAN  Data  Files:  None 

Requ i red  Subrout ine  s : None 

Method : Non-zero  linearly  interpolated  B-spline  coefficients,  N,  are 

computed  as  in  Equation  (4)  and  stored  in  the  array  CN.  A matrix  packing 
index,  LET,  is  computed  for  later  reconstruction  of  a full  vector. 

LFT  = S* FLOAT  (Ml -1)+1 , for  a non-periodic  basis.  For  a periodic  basis  the 
computation  of  LFT  is  made  modulo  (Ml-1). 


MBORD  (integer/ input) 


CS  (mixed/array/ input) 


K1  (integer/ input) 

I CLOSE  (integer/ input) 
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Remarks : 

(1)  The  combination  of  the  interpolation  algorithm  and  the  matrix 
packing  technique  requires  that  K1  * Ml. 

(2)  This  subroutine  is  used  in  the  same  manner  as  subroutine 
MKSPLN  and  all  the  "Remarks"  which  pertain  to  MKSPLN  also  apply  to  this 
subroutine.  In  general,  MBORD+1  non-zero  interpolated  coefficients  are 
produced  by  a call  to  PLSPLN,  rather  than  the  MBORD  non -zero  coefficients 
computed  by  a call  to  MKSPLN.  This  fact  must  be  accounted  for  in 
subsequent  calculations  and  prevents  the  direct  substitution  of  subroutine 
PLSPLN  for  MKSPLN. 
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FITTING  SUBROUTINES 
Subroutine  Name 


FITB01 

FITBS2* 


Page 
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FITB01  - Subroutine  Description 

Function:  To  compute  a least -squares -fit  transformation  matrix  for 

B-spline  curves  and  surfaces. 

Entry  Point:  FITB01 

Calling  Sequence:  CALL  FITB01(CS,CS,M1  ,MORIG,U,IOPT,IFS,SCR) 


Name  (Attributes) 

CS  (mixed/array/ input) 


Ml  (integer/ input) 

MORIG  (integer/ input) 

U (real/array/output) 

IOPT  (not  used) 

IFS  (integer /input) 


SCR  (array/scratch) 


COMON  Areas : None 

FORTRAN  Data  Files:  None 


Contents 

B-spline  coefficient  matrix  in 
packed  form.  Dimensioned 
CS(5, MORIG).  The  array  CS  is 
required  as  the  first  argument 
for  real  references  and  as  the 
second  argument  for  integer 
references . 

The  number  of  B-spline  basis 
functions  used  for  fitting. 

The  number  of  data  points  to  be 
fit. 

Least -squares  fit  transformation 
matrix.  Dimensioned  U (Ml , MORIG) . 

Flag  word.  If  zero,  CS  contains 
coefficients  for  a non-periodic 
or  open  B-spline  function.  If 
one,  CS  contains  coefficients  for 
a periodic  or  closed  B-spline 
function. 

A working  storage  area  KORE  words 
long,  where  K0RE=(M(Ml+l)*Ml)/2. 


Required  Subroutines : BCKSUB , BSMULT , DECOMP 
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Method : This  subroutine  computes  the  least-squares  transformation  [U] 

given  by  Equation  (8).  Subroutine  BSMULT  is  called  to  compute  the  product 
[S  ] [SI , where  the  coefficient  matrix  [S]  is  given  in  packed  form  in  array 
CS.  Subroutines  DECOMP  and  BCKSUB  are  called  to  solve  the  system  of 
equations  [S^S] [U1  = [S^]  for  [U] . 

Remarks : 

(1)  A program  stop  "STOP  67"  will  be  made  whenever  errors  occur 
T 

which  cause  the  matrix  [S  S]  to  be  singular.  This  should  occur  only  when 
CS,  Ml,  and  MORIG  are  not  compatible. 
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FITBS2  - Subroutine  Description 

Function:  To  determine  a polygon  or  two-dimensional  mesh  which  will 

induce  a B-spline  curve  or  surface  that  best  fits  a set  of  curve  or 
surface  data  points  in  the  least -squares  sense. 

Fjitry  Point:  FITBS2 

Calling  Sequence:  CALL  FITBS2(Q, P.MD.MORIG, NORIG, Ml, N1 ,ITYPE, RMS, 
5 IOPT,SCR,IFAIL) 


Calling  Arguments: 

Name  (Attributes) 

Q (real/array/input) 


P (real/array/output) 


MD, MORI G, NORIG  (integer/ 
input) 

Ml  (integer/ input/output) 


N1  (integer/ input/output) 


ITYPF.  (integer/ input) 


RMS, I OPT  (not  used) 
SCR  (array/scratch) 


Contents 

Rectangular  array  of  curve  or 
surface  data  points.  Dimensioned 
Q(MD,M0RIG, NORIG). 

Polygon  or  mesh  defining  a B- 
spline  curve  or  surface. 
Dimensioned  P(MD,M1,N1). 
Dimensions  of  array  Q. 

The  number  of  B-spline  basis 
functions  to  be  used  to  constiuct 
the  B-spline  curve  or  the  number 
of  basis  functions  to  be  used  in 
the  S-direction  of  a B-spline 
surface . 

The  number  of  B-spline  basis 
functions  to  be  used  in  the  t- 
direction  of  a B-spline  surface. 
Code  word  indicating  type  of  B- 
spline  curve  or  surface  (see 
Table  1). 

A working  storage  area  which  is 
KORL  words  long,  where 
K0RF.=M1  *MOR I G+MAXO ( K1  ,K2,K3) , 
K1=5*M0RIG+(M1*  (Ml  + 1) )/2 , 
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K2=Nl*NORIG+5*NORIG+ (Nl* (Nl+1 ) ) /2 , 
and 

K3=N1 *NORI G+MD*M1 *NORI G . 

IFAIL  (integer/output)  Flag  word.  A value  of  zero 

indicates  successful  completion. 

A value  of  1 indicates  that  the 
specified  combination  of  MORIG, 
NORIG,  Ml,  Nl,  and  I TYPE  is 
incompatible.  A value  of  2 
indicates  that  an  illegal  value 
of  ITYPF.  has  been  encountered. 

CONMON  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  BSCMAT , FITB01,  MULMAT 

Method : This  subroutine  will  evaluate  Equation  (7)  to  obtain  the 

polygon  P from  a set  of  curve  data  points  Q or  Equation  (11)  to  obtain  the 
mesh  [ P]  from  a set  of  surface  data  points.  Subroutine  BSCMAT  is  called 
to  compute  the  B-spline  coefficient  matrices  [S]  and  f T J . Next,  subroutine 
FITB01  is  called  to  compute  the  least -squares  transformations  [U]  and  [W] 
using  Equations  (8)  and  (12).  Then,  subroutine  MULMAT  is  called  to 
perform  the  final  multiplications  to  obtain  P,  or  [P] . 

Remarks : 

(1)  If  the  value  of  Ml  is  zero,  it  is  changed  to  the  value  of 
MORIG  and  the  transformation  [U]  is  replaced  by  the  identity  matrix. 
Similarly,  if  Nl  is  zero,  it  is  set  to  NORIG  and  [W]  is  replaced  by  the 
identity  matrix. 

(2)  The  arrays  Q and  P may  overlap  for  surface  fitting  applications 
whenever  Q is  not  required  for  subsequent  calculations.  The  arrays  Q and 

P may  not  overlap  for  curve  fitting  applications. 

(3)  The  transformation  matrices  [S]  and  [T] , as  computed  by- 
subroutine  BSCMAT,  are  based  on  a uniform  spacing  of  parameter  values 
between  zero  and  one.  A weighted  least -squares  fit  can  be  obtained  by 
providing  a substitute  for  BSCMAT  which  will  compute  a non -uni font 
transformation  matrix. 
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ABC  - Subroutine  Description 


Function:  To  determine  where  the  boundary  of  one  of  the  surfaces 

intersects  the  other  surface. 

Entry  Point:  ABC 

Calling  Sequence:  CALL  ABC(l  ,J ,K,ITOL,NGU') 

Calling  Arguments: 


Name  (Attributes) 
I (integer/ input) 


Contents 


J (integer/ input) 


K (integer/ input) 


Surface  whose  boundary  is  to  be 
used . 

Describes  which  part  of  boundary 
(of  surface  which  is  being 
considered  in  terms  of  its 
boundary)  is  to  be  used. 

1 = boundary  for  which  s is  constant. 

2 = boundary  for  which  t is  constant. 
Describes  which  part  of  boundary 
(of  surface  which  is  being  con- 
sidered in  terms  of  its  boundary) 
is  to  be  used. 

1 = parameter  which  is  to  be  constant 
will  be  made  0.0. 

2 = parameter  which  is  to  he  constant 
will  be  made  1.0. 

Control  parameter  for  iterative  non- 
linear minimization  subroutine. 

Number  of  attempts  to  be  made  to 
find  intersection  points  along  each 
edge  of  the  two  surfaces  during 
edge  analysis. 

boprxx,  stcon , bond,  sel 

FORTRAN  Data  Files:  None 


ITOL  (integer/ input) 
NOU  (integer/ input) 


COMON  Areas: 


Required  Subroutines:  OUTRNG , GINTR,  RF2SS,  RF2ST,  STR 

Method:  Use  iterative  solver. 
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CHKR  - Subroutine  Description 


Function:  To  evaluate  the  outcome  of  a call  to  the  iterative  solver 


Entry  Point : CHKR 


Calling  Sequence:  CALL  CHKR(I ,J,K,IRF,PARAM) 


Name  (Attributes) 
I (integer/ input) 


Contents 


Specifies  which  surface  is  of 
concern  for  solution  just 
attempted. 

0 = all  three  surfaces 

1 = Surface  i (i=l,2,  or  3) 
Specifies  which  type  of  conver- 
gence is  considered  acceptable 
from  solution  just  attempted. 

0 = Accept  EPS -convergence  (all 
residuals  less  than  the  specified 
tolerance)  or  NSIG -convergence 
(independent  variables  for  two 
consecutive  iterations  all 
agreed  to  within  the  specified 
tolerance) . 

1 = Accept  EPS -convergence  only. 
Tells  whether  solution  from 
iterative  solver  is  acceptable. 

0 = Pass;  1 = Fail. 

Termination  parameter  returned  by 
iterative  solver  which  relates 
what  happened  when  it  was  called. 
Parameters  (s^ ,tj»s, ,t^) . 
Dimensioned  PARAM(6) . 


J (integer/ input) 


K (integer/output) 


IRF  (integer/ input) 


PARAM  (real/array/input/ 
output) 

CONMON  Areas:  BOPRXX 

FORTRAN  Data  Files:  None 


Required  Subroutines:  OUTRNG 


EFF  - Function  Subroutine  Description 


Function:  To  build  termination  code  word  for  intersection  subroutines. 

Entry  Point:  EFF 

Calling  Sequence:  X = EFF(I,J) 


Calling  Arguments: 

Name  (Attributes) 

EFF  (real/ function/ output) 
I (integer/ input) 

J (integer/ input) 

CONMON  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  None 


Contents 

Termination  code  word. 

Number  used  to  indicate  where  in 
che  main  subroutine  a failure 
occurred . 

Number  which  relates  outcome  of 
a call  to  subroutine  GINTR. 
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GINTR  - Subroutine  Description 

Function:  To  iteratively  minimize  a set  of  nonlinear  functions. 

Entry  Point:  GINTR 

Calling  Sequence:  CALL  GINTR (RESFN ,FAR , ITOL ,PARVAL ,NRES ,NVAR , I FAIL) 


Name  (Attributes] 


Contents 


RESFN  (real/extemal  function/  Function  used  to  compute 


input] 

FAR  (real/input) 

ITOL  (integer/ input) 

PARVAL  (real/array/ input/ 
output) 


NRES  (integer/ input) 

NVAR  (integer/ input) 
IFAIL  (integer/output) 


residuals . 

Control  parameter  for  iterative 
solver. 

Control  parameter  for  iterative 
solver. 

The  independent  variables  for  the 
nonlinear  functions. 

As  input:  initial  guesses. 

As  output:  solution. 

Dimensioned:  PARVAL (NVAR) 

Number  of  dependent  variables 
(nonlinear  functions) . 

Number  of  independent  variables. 
Parameter  which  indicates  outcome 
of  a call  to  this  subroutine. 

-1  = iterative  solver  failed  to 
converge. 

0 = iterative  solver  converged 

with  all  function  values  less 
than  the  specified  tolerance. 

1 = iterative  solver  converged 

when  the  independent  variables 
for  two  consecutive  iterations 
all  agreed  to  within  a specified 
tolerance. 

7 = the  working  storage  required 

for  the  iterative  solver  was 
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not  available.  (Increase 
the  length  of  the  /BOPRXX/ 
COWON  block.) 


COMMON  Areas : RESCX , BOPRXX 

FORTRAN  Data  Files:  #6  (output/formatted) 

Required  Subroutines:  MKMARQ 

Method:  Determine  whether  sufficient  memory  is  available  and  call 

on  the  iterative  nonlinear  minimization  subroutine  MKMARQ. 


i 
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INT2S  - Subroutine  Description 

Function:  To  find  a sequence  of  points  on  the  curve  of  intersection 

of  two  surfaces. 

Entry  Point:  INT2S 

Calling  Sequence:  CALL  INT2S(IRESFN,FAR,IT0L,P0INTS,IMULT,IFAIL, 

NPOINT, ICOORD , KPLANE , CORVAL ) 

Calling  Arguments: 


Name  (Attributes) 
IRESFN  (unused) 
FAR  (real/ input) 


ITOL  (integer/ input) 


POINTS  (real/array/output) 


IMJLT  (unused) 

IFAIL  (integer /output) 


■ITOL 


Contents 

Control  parameter  for  iterative 
nonlinear  minimization  subroutine; 
selects  how  two  methods  used  are 
to  be  mixed.  In  absence  of  any 
special  preference,  should  be  set 
at  1.0. 

Control  parameter  for  iterative 
nonlinear  minimization  subroutine. 
A solution  is  considered  to  have 
been  reached  if  either  of  the 
following  occurs: 

(1)  all  functional  values  <10 

(2)  independent  variables  all 
agree  to  ITOL  significant  figures 
in  any  two  consecutive  iterations. 
Polygon  of  points  on  the  curve  of 
intersection.  Dimensioned 
POINTS (3, NPOINT) . 

Value  Meaning 

0 Curve  of  intersection 

successfully  found. 
lOOOi+j  No  curve  of  intersection 
has  been  found;  see 
"Remarks". 
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curve  of  intersection. 

ICOORD  (integer/ input)  Value  Meaning 

0 Do  not  substitute  a coordinate 
plane  for  the  second  surface. 

1 Do  substitute  a coordinate 
plane  for  the  second  surface. 

KPLANF.  (integer/ input)  Descriptor  for  coordinate  plane 

being  substituted  for  one  of  the 
surfaces;  used  only  if  IC00RD=1. 

1 = Plane  of  constant  x 

2 = Plane  cf  constant  y 

3 = Plane  of  constant  z 

CORVAL  (real/input)  Descriptor  for  coordinate  plane 

being  substituted  for  one  of  the 
surfaces;  used  only  if  IC00RD=1; 
value  of  constant  coordinate. 

CONMON  Areas:  BOND,  BOPRXX,  RESCX,  SEL,  STCON,  TIE 

FORTRAN  Data  Files:  None 

Required  Subroutine:  ABC,  EFF,  GINTR,  OUTRNG,  RF2SS,  RF2ST 

Method:  Iterative  nonlinear  minimization  scheme. 

Remarks : 

(1)  Regarding  IFAIL:  If  i is  1,  2,  or  3,  failure  occurred  during 

attempt  to  lay  out  a curve  of  intersection  between  two  points.  If  i = 1: 
failed  to  converge;  if  i = 2,  3:  point  out  of  range  on  one  of  the  surfaces; 
if  i = 4:  problem  is  inadmissible,  j is  always  1 for  i = 2,  3,  or  4. 

If  i = 1,  j indicates  outcome  from  iterative  solver  as  follows:  If  j = 1: 
iterative  solution  converged  with  all  residuals  less  than  a specified 
tolerance;  if  j = 2:  iterative  solution  converged  to  a stationary  point; 
if  j = 3:  iterative  solution  failed  to  converge;  and  if  j = 4:  insufficient 
memory  for  iterative  solver. 

(2)  Regarding  the  following  CGNMON  block:  CCMDN/ BOPRXX/ LBOPR 1 , 

LB0PR2 , I.BIJFF,  JBUFF,  BUFF(l),  MBOPR(8,9),  BUFFX(87) 

The  information  described  below  must  be  entered  into  this 
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block  before  INT2S  is  called.  The  main  function  of  this  block  is  to 
provide  a storage  area  for  the  definitions  of  the  two  surfaces.  The  main 
subdivisions  of  this  block  are  described  in  Table  2.  Array  MBOPR  contains 
data  about  the  two  surfaces.  The  surface  data  region  is  described  in 
Table  3.  Array  BUFF  must  contain  the  grid  of  defining  points  for  each  of 
the  two  surfaces.  These  grids  may  be  entered  into  BUFF  consecutively. 

BlfFF  is  also  a working  storage  area. 

TABLE  2 - THE  MAIN  SUBDIVISIONS  OF  COMON  BLOCK  BOPRXX 


Name 

LB0PR1 , LB0PR2 

LBUFF 

JBUFF 

BUFF 

MBOPR 


BUFFX 


Description 

Dimensions  of  MBOPR. 

Length  of  BUFF. 

Next  available  location  in  BUFF. 

Dummy  working  storage  area  to  permit  both 
MBOPR  and  BUFFX  to  be  assigned  dynamically. 
Data  index  containing:  surface  data 
pointers,  spatial  guess*,  and  transfor- 
mation matrix  data  pointers*. 

Working  storage  area;  includes  surface 
defining  points  and  transformation 
matrices*. 


This  information  is  not  required  for  subroutine  INT2S. 


TABLE  3 - SURFACE  DATA  AREA  OF  MBOPR 


Location* 

Descript  ion 

MB0PR(2, I) 

Surface  type;  see  Table  1 

MBOPR (3,1) 

Pointer  to  first  word  of  sur- 
face definition  in  BUFF. 

MB0PR(4,I) ,MROPR(5,T) 

Dimensions  of  defining  mesh 
for  surface. 

* I is  1,  2,  or  3 and  denotes  which  surface  is  being  described. 


•<- 
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INT2SX  - Subroutine  Description 


Function:  To  calculate  the  location  of  points  on  the  curve  of  inter- 

section of  two  B-spline  surfaces. 

Entry  Point:  INT2SX 

Calling  Sequence:  CALL  INT2SX(P1jM1 ,N1 , ITYP1 ,P2 ,M2 ,N2 , ITYP2 , 

ICOORD , KPLW^-ORVAL  .POINTS  .NPOINT  , ITOL , 
IMOVE.IFAIL) 


Calling  Arguments: 

Name  (Attributes) 

Pi  (real/array/input) 

Mi,Ni  (integer/ input) 
ITYPi  (integer/ input) 

ICOORD  (integer/ input) 


KPLANE  (integer/ input) 


CORVAL  (real /input) 

POINTS  (real/array/output) 


Contents 

Mesh  defining  the  i^  B-spline 
surface.  Dimensioned  Pi(3,Mi,Ni). 
Dimensions  of  array  Pi. 

Code  word  indicating  the  type  of 
B-spline  surface  stored  in  Pi 
(see  Table  1) . 

Flag  word.  If  IC00RD=1,  inter- 
section problem  will  be  solved 
using  a specified  coordinate 
plane  for  the  second  surface.  If 
ICOORD=0,  an  arbitrary  surface  is 
assumed  to  be  present  as  the 
second  surface . 

Code  word  indicating  orientation 
of  coordinate  plane,  if  used. 

If  KPLANE=1,  an  X-plane  will  be 
used. 

If  KPL\NE=2,  a Y-plane  will  be 
used . 

If  KPLANE=3,  a Z -plane  will  be 
used. 

Coordinate  value  for  plane  used. 
Array  of  points  calculated  on  the 
curve  of  intersection. 

Dimensioned  POINTS (3, NPOINT) . 
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NPOINT  (integer/ input) 


ITOL  (integer/ input) 


IMOVF.  (integer/ input) 


COMMON  Areas : 


BOPRXX 


Number  of  points  to  be  calculated 
on  the  curve  of  intersection. 

Must  be  - 2. 

Convergence  tolerance  for 
iterative  solver.  Procedure  will 
terminate  when  an  additional 
iteration  docs  not  change  the  ITOL 
most  significant  decimal  digits 
of  the  solution. 

Flag  word.  If  IM0VE=1 , surface 
meshes  will  be  transferred  to 
COMMON  block  /BOPRXX/  prior  to 
solution.  If  IM0VF=O,  only  the 
memory  addresses  of  the  surface 
meshes  will  be  stored  in  CONMON 
block  /BOPRXX/. 

Flag  word.  A value  of  IFAIL 
which  is  less  than  1000  indicates 
a successful  completion  of  the 
procedure.  A value  which  is 
1000  or  greater  indicates  a 
failure  condition  (see  description 
of  subroutine  INT2S  for  details). 


Function 

Communication  with  subroutine  INT2S  (see  description 
of  subroutine  INT2S  for  structure) . 


IFAIL  (integer/output) 


FORTRAN  Data  Files:  None 

Required  Subroutines:  INT2S,  LOCF 

Method : This  subroutine  is  a machine  dependent  driver  for  solving  the 

two-surface  intersection  problem  outside  of  the  G-PRIME  environment. 
Operands  are  moved  to  COfWON  block  /BOPRXX/  and  then  subroutine  INT2S  is 
called  to  calculate  the  solution. 
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Remarks'- 

(1)  Subroutine  LOCF,  used  to  obtain  the  memory  addresses  of  the 
surface  operands,  is  machine  dependent.  Use  of  this  subroutine  may  be 
avoided  by  setting  IMOVE  to  1 and  increasing  the  length  of  COMMON  block 
/BOPRXX/  to  accommodate  copies  of  the  surface  operands. 


INT3S  - Subroutine  Description 

Function:  To  find  a point  of  intersection  of  three  surfaces. 


INT3S (IRESFN , FAR , ITOL , PAKVAL , XYZ , INULT , I FAIL) 


Fjitry  Point:  INT3S 

Calling  Sequence:  CALL 
Calling  Arguments: 

Name  (Attributes) 
IRESFN  (integer/ input) 


FAR  (real /input) 

ITOL  (integer/ input) 


PARVAL  (real/array/ input/ 
output) 


Contents 

Flag  word.  If  IRESFN=1  or 

IRESFN=2  solve,  starting  with 

parameter  guesses,  using  residual 

functions  RF3S01  or  RF3S02, 

respectively.  If  IRESFN=S  or 

IRESFN'=6  locate,  then  solve  using 

residual  functions  RF3S01  or 

RF3S02,  respectively. 

Control  parameter  for  iterative 

nonlinear  minimization  subroutine. 

Control  parameter  for  iterative 

nonlinear  minimization  subroutine 

used  for  "solve".  A solution  is 

considered  to  have  been  reached  if 

either  of  the  following  occurs: 

- TTOI 

(1)  all  functional  values  • 10  J 

(2)  independent  variables  all 
agree  to  ITOL  significant  figures 
in  any  two  consecutive  iterations. 
Independent  variable  for  nonlinear 
solver.  Dimensioned  PARYAL(6) . 
Input:  initial  guesses  for  the 
parameters  s^  ,t^  ,s., ,t? ,s. , and  t.. 
If  doing  "Locate -Solve",  use  0.S 
for  all  six  values.  If  doing 
"Solve",  use  best  estimate  of 
parameter  values  at  point  of 
intersect  ion. 


49 


XYZ  (real/array/output) 


IMULT  (integer/ input) 


IFAIL  (integer/output) 


Output:  if  a solution  has  been 
found,  PARV/M,  will  contain  the 
parameter  values  s^  ,t j ,s,  ,t-,  ,s^ , 
and  t,,  associated  with  the  point 
of  intersection  of  the  three 
surfaces . 

Coordinates  of  point  of  inter- 
section of  the  three  surfaces, 
if  found.  Dimensioned  XYZ(3) . 
Parameter  which  determines 
whether  or  not  the  program  is  to 
make  any  additional  effort  at 
finding  a solution  when  the  first 
try  fails.  fSec  "Remarks".')  If 
IMIJLT=1,  make  additional  effort. 

If  IMIJLT > 1 , do  not  make  additional 
effort . 

If  IFAIL=0,  an  intersection  point 
has  been  found.  If  IFAIL=1000i+j , 
an  intersection  point  has  not 
been  found;  see  "Remarks". 


COMO.M  Areas : ARB,  ROPRXX,  MULT,  RESCX 

FORTRAN'  Data  Files:  None 

Required  Subroutines:  CHKR,  F.FF,  GINTR,  MF.K,  RF3S01 , RF3S02,  XLiXIl, 

XL0CI2,  XI.0CI3 

Method:  Iterative  nonlinear  minimization  scheme. 

Remarks : 

(1)  If  initial  guess  is  spatial  (x,y,z),  it  must  be  stored  in 
MBOPR(3,4),  MBOPR(4, 4) , MR0PR(S,4)  in  CCMION  block  /ROPRXX/  before  calling 
the  subroutine.  If  initial  guess  is  parametric  (Sj  ,t  { .s.,  ,t,  ,s,  ,t ^ . it 
must  be  put  in  an  array  which  will  become  PARVAL(l)  through  PAR\ AL(6)  when 
this  subroutine  is  called. 

(2)  Regarding  IMULT:  When  IMULT=1,  the  "additional  effort"  made 
for  spatial  initial  guesses  and  for  parametric  initial  guesses  is  quite 
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different.  When  the  initial  guess  is  parametric,  the  "additional  effort" 
consists  of  substituting  various  sets  of  parametric  initial  guesses 
(s^  ,t^  ,St ,t? .Sj,t j)  for  the  parametric  initial  guess  supplied  by  the  user. 
When  the  initial  guess  is  spatial,  the  spatial  initial  guess  supplied  by 
the  user  remains  the  same;  the  "additional  effort"  is  made  during  the 
"locate"  process  for  one  of  the  three  surfaces  and  consists  of  trying  a 
new  initial  guess  for  the  parameters  (s^,t^)  for  the  surface  on  which  the 
"locate"  is  being  attempted.  No  "additional  effort"  is  made  during 
the  "solve"  phase. 

(3)  Regarding  IFAIL:  Generally,  i indicates  the  place  in  INT3S 

at  which  failure  occurred;  j indicates  the  outcome  of  the  most  recent 
call  to  GINTR  (which  uses  the  iterative  nonlinear  minimization  process). 

i_  Meaning 

1 Failed  during  "solve" 

2 Failed  during  "locate"  for  first  surface 

3 Failed  during  "locate"  for  second  surface 

4 Failed  during  "locate"  for  third  surface 

5 Failed  during  "solve"  of  "locate-solve" 

Meaning 

1 Iterative  solution  converged  with  all  residuals 

-TTDI 

less  than  10  '.  Usually  indicates  that  the 

solution  point  is  outside  standard  surface 
definition  limits. 

2 Iterative  solution  converged  to  a stationary  point. 

3 Iterative  solution  failed  to  converge. 

4 Insufficient  memory  for  iterative  solver. 

(4)  Regarding  the  following  C(*WON  block:  CONMON/BOPRXX/LBOPR1 , 
LBOPR2 , LBUFF , JBUFF , BUFF ( 1 ) ,MBOPR(8,9)  ,RIJFFX(87) 

The  information  described  below  must  be  entered  into  this 
block  before  INT3S  is  called.  The  main  purpose  of  this  block  is  to  pro- 
vide a storage  area  to  contain  the  definitions  of  the  three  surfaces.  The 
main  subdivisions  of  this  block  are  described  in  Table  2.  Array  MROPR 
contains  three  main  storage  areas: 

(1)  data  about  the  three  surfaces 


51 


(2)  a spatial  guess  for  the  intersection  point 

(3)  data  about  the  transformation  matrices  (which  are 
required  whenever  the  piecewise -linear  approximation  for  the  three 
surfaces  is  used) . 

The  surface  data  area  is  described  in  Table  3.  The  spatial  guess  must  be 
put  into  MBOPR(3,4),  MBOPR(4,4),  and  MROPR(5 ,4) ; however,  the  spatial 
guess  must  be  entered  into  these  locations  as  real  numbers.  The  trans- 
formation matrix  data  area  is  described  in  Table  4. 

TABLE  4 - TRANSFORMATION  MATRIX  DATA  AREA  OF  MBOPR 


Location* 

Description 

MBOPR(3 , I) 

Pointer  to  first  word  of  trans- 

formation  matrix  in  BUFF. 

MBOPR (4,1) 

Number  of  defining  points. 

MROPR (5,1) 

Number  of  output  points. 

* I is  8 = open,  9 = closed. 


Array  BUFF  must  contain: 

(1)  the  grid  of  defining  points  for  each  of  the  three 
surfaces,  and 

(2)  the  transformation  matrices. 

The  transformation  matrices  may  be  generated  by  calling  subroutine  BSCMAT 
as  follows: 


CALL  BSCMAT  ( BUFF  ( I ) , . J , K . L) 

Pointer  to  first  word  of  transformation  matrix  in  BUFF,  i.e., 
the  location  in  BUFF  where  the  user  wants  the  transformation 
matrix  to  begin. 

Number  of  defining  points 
Number  of  output  points 
0 = Open;  1 = Closed 
Blocks  of  data  should  be  entered  into  BUFF  consecutively;  .TBUFF  should  be 
updated  each  time  a new  block  of  data  is  entered  into  BUFF.  INT3S  requires 
87  words  of  working  storage  in  the  BUFFX  array. 


where  l 


J 

K 

L 
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INT3SX  - Subroutine  Description 
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Function:  To  calculate  the  location  of  a point  of  intersection  of 

three  B-spline  surfaces. 

Entry  Point:  INT3SX 

Calling  Sequence:  CALL  INT3SX(P1 ,M1  ,N1 , ITYP1 ,P2 ,M2 ,N2 , ITYP2 ,P3 , 

M3, N3,ITYP3,SPOINT, PARVAL, XYZ, ITOL, LOCATE, 
IMOVE , I FAIL) 


Calling  Arguments: 

Name  (Attributes) 

Pi  (real/array/input) 

Mi,Ni  (integer/ input) 
ITYPi  (integer/ input) 


SPOINT  (real/array/ input) 


PARVAL  (real/array/input/ 
output) 


XYZ  (real/array/output) 
ITOL  (integer/ input) 


Contents 

Mesh  defining  the  i^h  B-spline 
surface.  Dimensioned  Pi(3,Mi,Ni). 
Dimensons  of  array  Pi. 

Code  word  indicating  the  type  of 
B-spline  stored  in  Pi  (see  Table 
1). 

Coordinates  of  spatial  point  close 
to  the  desired  intersection  point. 
Dimensioned  SP0I\T(3). 

On  input,  parametric  guesses  of 
the  location  of  the  desired  point 
on  each  of  the  surfaces.  If 
unknown,  values  of  .5  are  a 
reasonable  choice.  On  output, 
contains  the  parametric  location 
of  the  intersection  for  each  of  the 
three  surfaces.  Dimensioned 
PARVAL (6) . 

Spatial  coordinates  of  the  inter- 
section point.  Dimensioned  XYZ(3). 
Convergence  tolerance  for  iterative 
solver.  Procedure  will  terminate 
when  an  additional  iteration  does 
not  change  the  ITOL  most  significant 
decimal  digits  of  the  solution. 
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LOCATE  (integer/ input) 


Flag  word.  If  liOCATE=l, 
coordinates  in  array  SPOINT  will 
be  used  to  determine  a new  set  of 
initial  parameter  guesses  for 
iterative  solving  procedure.  If 
LOCATE=0,  parameter  values  in 
PARVAL  will  be  used  as  ini t i ° 1 
guesses . 

WOVE  (integer/ input)  Flag  word.  If  IM3VE=1,  surface 

meshes  will  be  transferred  to 
COWON  block  /BOPRXX/  prior  to 
solution.  If  IMOVE=0,  only  the 
memory  addresses  of  the  surface 
meshes  will  be  stored  in  COMON 
block  /BOPRXX/. 

I FAIL  (integer /output)  Flag  word.  A zero  value  of  IFAIL 

indicates  a successful  completion 
of  the  procedure.  All  other 
values  indicate  a failure  condition 
(see  description  of  subroutine 
INT3S  for  details). 

COMON  Areas: 

Name  Function 

BOPRXX  Communication  with  subroutine  INT3S  (see 

description  of  subroutine  INT3S  for  structure). 

FORTRAN  Data  Files:  None 

Required  Subroutines:  INT3S,  LOCF 

Method:  This  subroutine  is  a machine  dependent  driver  for  solving  the 

three-surface  intersection  problem  outside  of  the  G-PRIME  environment. 
Operands  are  moved  to  COWON  block  /BOPRXX/  and  then  subroutine  INT3S  is 
called  to  calculate  the  solution. 

Remarks : 

(1)  Subroutine  LOCF,  used  to  obtain  the  memory  addresses  of  the 
surface  operands,  is  machine  dependent.  Use  of  this  subroutine  may  be 
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avoided  by  setting  IMOVE  to  1 and  increasing  the  length  of  COMON 
block  /BOPRXX/  to  accommodate  copies  of  the  surface  operands. 


MEK  - Subroutine  Description 

Function:  Select  a new  parameter  guess  for  iterative  solver. 

Entry  Point:  MEK 

Calling  Sequence:  CALL  MEK(I ,J,K,L, PARAM) 

Calling  Arguments: 

Name  (Attributes) 

I (integer/ input/output) 


Contents 


J (integer/ input) 
K (integer/ input) 


L (integer/output) 


PARAM  (real/array/input/ 
output) 

COfMON  Areas:  ARB,  MULT 

None 
None 


Indicates  which  parameter  guess 
is  currently  being  used. 

Maximum  number  of  parameter 
guesses  allowed. 

Indicates  which  surfaces  are 
involved. 

0 = all  three  surfaces 

1 = Surface  i (i  = 1,  2,  or  3) 

0 = A new  parameter  guess  was 
made;  this  implies  that  another 
solution -attempt  may  be  made. 

1 = No  more  parameter  guesses  are 
available;  a new  parameter  guess 
was  not  made. 

Parameters  (s^ ,t^ ^^.s^.t-j) . 
Dimens i oned  PARAM (2,3). 


FORTRAN  Data  Files: 


Required  Subroutines: 

Method:  Select  additional  parameter  guesses  from  table. 
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OUTRNG  - Subroutine  Description 


Function:  To  process  specific  values  for  the  parameters  (s,t)  for  a 


surface  as  follows:  (a)  if  parameter  is  "open",  determine  whether  its 
value  is  "out  of  range";  (b)  if  parameter  is  "closed",  bring  its  value 
back  into  range. 

F.ntry  Point : OUTRNG 

Calling  Sequence:  X = OUTRNG ( ITYPE , PARPAJ^, EPS) 

Calling  Arguments: 


Name  (Attributes) 

OUTRNG  (logical/output) 


Contents 


ITYPE  (integer/ input) 
PARPAR  (real/array/input/ 
output) 

EPS  (real /input) 


Status  indicator.  If  .TRUE.,  at 
least  one  parameter  value  is  out 
of  range.  If  .FALSE.,  all 
parameter  values  are  in  the  range 
[0,1]. 

Entity  type;  see  Table  1. 

Surface  parameters  (s,t). 
Dimensioned  PARPAR(2). 

Tolerance  value  used  in  "out-of- 
range"  determination. 


COMON  Areas:  None 

FORTRAN  Data  Files:  None 


Required  Subroutines:  None 

Method:  Comparison  using  specified  tolerance.  , 

Remarks : 

(1)  When  a parameter  is  "out  of  range",  it  lies  outside  the 


interval  [0,1].  The  tolerance  value,  EPS,  is  used  to  expand  the  interval 
by  a small  amount. 


57 


RF2SS  - Function  Subroutine  Description 


Function:  To  compute  the  residual  functions  needed  to  solve  for  one 

point  on  a curve  of  intersection  of  two  surfaces. 

Entrv  Point:  RF2SS 


Calling  Sequence:  Z = RF2SS(X,I) 

Calling  Arguments: 

Name  (Attributes) 

RF2SS  (real /function/output) 

X (real/array/input) 


Contents 


I (integer/ input) 


Value 

1 

2 

3 


Residual  value. 

Parameters  for  the  two  surfaces; 
dimensioned  X(3);  X (1)  and  X(2) 
are  the  parameters  (s  and  t)  for 
the  surface  being  considered  in 
its  entirety;  X(3)  is  the 
variable  parameter  (t)  for  the 
other  surface. 

Indicates  which  of  the  three 
residual  functions  is  to  be 
computed . 

Meaning  (See  Remark  3) 

Xi(Si.ti)  ■ x2(s0>t2) 

yi(srV  " y2*'S0’t2'* 

Mvv  * z2^S0,t2-' 


COMMON  Areas:  SEL,  STCON,  RESCX,  BOPRXX 
FORTRAN  Data  Files:  None 

Required  Subroutines:  BSEVL1 

Method:  Compute  residual  functions  from  evaluation  of  spatial 

coordinates . 

Remarks : 

(1)  Regarding  COWON  block  STCON:  COWON/ STCON/ SCON , TCON 

SCON  must  contain  s^. 

(2)  Regarding  COWON  block  SEL:  COWON/SEL/KA , KB 

KA  denotes  the  surface  (Surface  1 or  Surface  2)  which  is 
to  be  considered  in  its  entirety. 
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KB  denotes  the  surface  (Surface  1 or  Surface  2)  which  is 
to  be  considered  only  in  terms  of  a curve  of  constant  s. 

(3)  s^  and  t^  are  the  parameters  for  the  surface  which  is  being 
considered  in  its  entirety. 

Sq  (a  constant)  and  t0  are  the  parameters  for  the  other 

surface. 


The  spatial  coordinates  for  the  surface  which  is  being 
considered  in  its  entirety  are  given  by 

^1  ^1  ’^1^ 
y1Cs1,t1) 

^i  ^i  ’^i^ 


The  spatial  coordinates  for  the  other  surface  are  given  by 
x2(so’t2^ 
y 2 


1 


f I 
i 


t » 
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RF2ST  - Function  Subroutine  Description 


To  compute  the  residual  functions  needed  to  solve  for  one 
of  intersection  of  two  surfaces. 


Z = RF2ST(X, I) 


Function: 
point  on  a cun/e 

Entry  Point:  RF2ST 

Calling  Sequence: 

Calling  Arguments: 

Name  (Attributes) 

RF2ST  (real/function/output) 
X (real /array/ input) 


I (integer/ input) 

Value 
1 

2 

3 


Contents 
Residual  value. 

Parameters  for  the  two  surfaces; 
dimensioned  X(3);  X(l)  and  X(2) 
are  the  parameters  (s  and  t)  for 
the  surface  which  is  being 
considered  in  its  entirety;  X ( 3) 
is  the  variable  parameter  (s) 
for  the  other  surface. 

Indicates  which  of  the  three 
residual  functions  is  to  be 
computed . 

Meaning  (See  Remark  3) 

Xl^l’V  ” (s2 

yj(s^,t|)  * ^2^s2 ,t:0^ 

2.(S|,t^)  * *1q) 


C0M10N  Areas:  SEE,  STCON,  RESCX , BOPRXX 
FORTRAN  Data  Files:  None 

Required  Subroutines:  BSEYL1 

Method:  Compute  residual  functions  from  evaluation  of  spatial 

coordinates. 

Remarks : 

(1)  Regarding  C.ONMON  block  STCON:  CONWON/STCON/SCON , TCON 

TCON  must  contain  t„. 

(2)  Regarding  CONMON  block  SEL:  COMMON /SEL/KA,  KB 

KA  denotes  the  surface  (Surface  1 or  Surface  2)  which  is  to 
be  considered  in  its  entirety. 
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KB  denotes  the  surface  (Surface  1 or  Surface  2)  which  is 
to  be  considered  only  in  terms  of  a curve  of  constant  t. 

(3)  s^  and  t^  are  the  parameters  for  the  surface  which  is  being 
considered  in  its  entirety. 

s^  and  tp  (a  constant)  are  the  parameters  for  the  other 

surface. 

The  spatial  coordinates  for  the  surface  which  is  being 
considered  in  its  entirety  are  given  by 

^1 (s i > t^ ) 

y1(s1,t1) 

^i (Si , tj) 

The  spatial  coordinates  for  the  other  surface  are  given  by 
x2(s2.t0) 

y2(s2’V 

Z2(‘S2’t0^ 


h" 
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RF3S01  -‘Function  Subroutine  Description 


. 


Function:  To  compute  the  spatial  coordinates  on  three  surfaces 

corresponding  to  a given  set  of  parameters  and  then  to  compute  the 
residual  functions  needed  to  solve  for  a point  of  intersection. 


Fjitrv  Point : RF3S01 

Calling  Sequence:  Z = RF3S01(X,I) 


Calling  Arguments: 

Name  (Attributes) 

RF3S01  (real/function/ 
output) 

X (real/array/input) 

I (integer/input) 

Value 

1 

2 

3 

4 

5 

6 

where  the  spatial  coordinates  of  a 


Contents 
Residual  value. 

Parametric  coordinates  (s^,t^,s9, 
t7,s,,t3)  for  the  three  surfaces. 
Dimensioned  X(6) . 

Indicates  which  of  the  six  resi- 
dual functions  is  to  be  computed. 
Meaning 

xl^sl ’t^  ~ x9(s2>t-,) 

yi(sl’V  ' y2(s2,1:2^ 

Zjfs^tj)  - z2(s2,t2) 
x2(s2,t2)  - x3(s3,t3) 

y2^S2,t:2^  " y3^S3’t3^ 
z2(s2  »t?)  - “3(,s3»i^) 
point  on  Surface  i are  given  by 


COMMON  Areas:  RFSCX , BOPRXX 

FORTRAN  Data  Files:  None 

Required  Subroutines:  BSFVL3 

Method:  Compute  residual  functions  from  evaluation  of  spatial 

coordinates. 
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RF3S02  - Subroutine  Description 


Function:  To  compute  the  spatial  coordinates  on  three  surfaces 

corresponding  to  a given  set  of  parameters  and  then  to  compute  the 
residual  functions  needed  to  solve  for  a point  of  intersection. 

FJitry  Point:  RF3S02 

Calling  Sequence:  Z = RF3S02(X,I) 

Calling  Arguments: 

Name  (Attributes)  Contents 

RF3S02  (real/function/output)  Residual  value 

X (real/array/ input)  Parametric  coordinates  (s^,t^,s 

t_,,s^,t^)  for  the  three  surface 
Dimensioned  X(6). 

I (integer/ input)  Indicates  which  of  the  six 

residual  functions  is  to  be 
computed. 


Value 

Meaning 

1 

X1^S1>t1) 

- x7(s7,t2) 

2 

y j j 

y2(s2'tJ 

3 

zi.(sr  V 

-7(s7,t7) 

4 

x2^s2,t:2^ 

X 2 (s  J 

5 

^2  ^s2  * ^2^ 

■ >'3(s3't3) 

6 

z2fs2,t2') 

- z3(s3,t3) 

where  the  spatial  coordinates  of  a point  on  Surface  i are  given  by 
XjfSj.tj),  yi(si,ti),  and  zi(s.,ti). 

COMMON  Areas:  RESCX , BOPRXX 

FORTRAN  Data  Files:  None 

Required  Subroutines:  BSEVL1 

Method : Compute  residual  functions  from  evaluation  of  spatial 

coordinates . 
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STR  - Subroutine  Description 

Function:  To  build  a file  of  points  for  which  (1)  the  boundary  of 

Surface  1 intersects  Surface  2,  or  (2)  the  boundary  of  Surface  2 inter- 
sects Surface  1. 

Entry  Point:  STR 

Calling  Sequence:  CALL  STR(BPP,EPSIL) 

Calling  Arguments: 

Name  (Attributes)  Contents 

BPP  (real/array/ input)  Point  for  which  (1)  the  boundary 

of  Surface  1 intersects  Surface  2, 
or  (2)  the  boundary-  of  Surface  2 
intersects  Surface  1;  given  in 
terms  of  parameters  (s.  ,tj  .s-,  ,t,) . 
Dimensioned  BPP(2,2). 

EPSIL  (real/input)  Tolerance  value  used  for  checking 

whether  a point  is  already  in  the 
file;  if  all  four  parameters 
agree  to  within  the  specified 
tolerance  value  with  the  parameters 
for  a point  already  in  the  file, 
no  new  point  is  entered  in  the  file. 

COMMON  Areas:  BOPRXX,  TIE,  BOND 

FORTRAN  Data  Files:  None 

Required  Subroutines:  None 

Method:  Compare  the  newly  found  point  to  the  points  already  in  the 

file.  If  the  new  point  is  different  from  all  the  points  already  in  the 
file,  enter  the  new  point  in  the  file. 

Remarks : 

(1)  Regarding  COWON  block  TIE:  COWON/TIE/AB(2 ,2 ,8) 

Array  AB  contains  the  aforementioned  file  of  points. 
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XLOCIl  - Function  Subroutine  Description 

Function:  To  compute  the  residual  functions  needed  to  perform  the 

"locate"  process  (i.e.,  finding  the  point  closest  to  the  spatial  gue  s) . 
Fntry  Point:  XLOCI 1 

Calling  Sequence:  Z = XL0CI1(X,I) 


Name  (Attributes)  Contents 

XLOCIl  (real/function/output)  Residual  value 

X (real /array/ input)  Parameters  (s,t)  for  the  surface 

on  which  the  "locate"  is  being 
attempted.  Dimensioned  X(2). 

I  (integer/ input)  Coordinate  direction  currently 

being  called  for:  1 = x 

2 = y 

3 = z 

CQNttON  Areas:  BOPRXX , RESCX 

FORTRAN  Data  Files:  None 

Required  Subroutines:  BSEVL1 

Method:  Compute  difference  (in  coordinate  direction  specified) 

between  point  on  surface  and  spatial  guess. 

Remarks : 

(1)  This  function  has  three  entry  points.  When  the  "locate" 
process  is  to  be  performed  for  Surface  i,  entry  point  XLOCI i should  be 
used  (i  is  1 , 2,  or  3) . 
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UTILITY  subroutines 
Sub  rout  ine  Name 


Page 


BCKSUB 

67 

BSMULT 

68 

DECOMP 

69 

MKMARQ 

71 

MULMAT 

74 

NOBSE3 
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BCKSUB  - Subroutine  Description 

Function:  To  perform  a vector  forward  backward  substitution  using  a 

Cholesky  factor  matrix. 

Putty  Point : BCKSUB 

Calling  Sequence:  CALL  BCKSUB (S,X,B,N) 

Calling  Arguments: 

Name  (Attributes)  Contents 

S (real/array/ input)  Lower  triangular  factor  matrix, 

stored:  S(1 ,1) ,S(2,1) , . . . , 

S(N ,N-1) ,S(N ,N) . Length  is 
(N(N+l)/2). 

X (real/  array/output)  Solution  vector.  Length  N. 

B (real /array/ input)  Right-hand  side  vector.  Length  N. 

N (integer/ input)  Order  of  matrix  S. 

QMDN  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  None 


T 

Method:  The  equation  [S][S  ]X  = B is  solved  for  X by  first  solving 

[S]Y  = B for  Y and  then  solving  [ST]X  = Y for  X,  where  [S]  is  a lower 
triangular  Cholesky  factor  matrix. 

Remarks : 

(1)  This  subroutine  is  used  with  subroutine  DFCOMP  to  solve  the 
equation  [ A ] X = B.  See  description  of  subroutine  DECOMP  for  details. 

(2)  The  arrays  B and  X may  overlap  whenever  B is  not  required  for 
subsequent  calculations. 
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BSMULT  - Subroutine  Description 


j 


- 


f 

t 


T 

To  compute  the  matrix  product  [S  ][S],  storing  only  the 

[S]  must  be  a B-spline 

packed  fom. 


C\LL  BSMULT (CS ,CS ,M1  , MORIG,  IFS, STSLTR) 


Function: 

lower  triangular  portion  of  the  result, 
coefficient  matrix  stored  in 
F.ntrv  Point:  BSMULT 

I -d. — . — 

Calling  Sequence: 

Calling  Arguments: 

Name  (Attributes) 

CS  (mixed/array/ input) 


Ml  (integer/ input) 
MORIG  (integer/ input) 
IFS  ( integer/ input) 


STSLTR  (real/array/output) 


Contents 

B-spline  coefficient  matrix  in 
packed  form.  Dimensioned  CS(5, MORIG). 
The  array  CS  is  required  as  the 
first  argument  for  real  references 
and  as  the  second  argument  for 
integer  references. 

Row  length  of  full  [SJ  matrix. 

Column  length  of  full  [S]  matrix. 

Flag  word.  If  zero,  CS  contains 
coefficients  for  a non-periodic 
or  open  B-spline  function.  If 
one,  CS  contains  coefficients  for 
a periodic  or  closed  B-spline 
function. 

The  product  matrix  stored  in 
lower  triangular  form. 


COMMON  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  None 

Method : Fach  row  of  CS  is  unpacked  and  its  contributions  to  the 

product  are  accumulated  in  the  array  STSLTR. 

Remarks : 


(1)  The  form  of  the  product  array  STSLTR  is  compatible  with  the 
requirements  of  the  Cholesky  decomposition  subroutine,  DFCOIP. 


I 
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DECCMP  - Subroutine  Description 

Function:  To  perform  Cholesky  decomposition  of  a real,  symmetr  c 

matrix  when  diagonal  and  lower  triangle  of  matrix  are  given. 

Entry  Point:  DECCMP 

Calling  Sequence:  CALL  DECCMP (A, S, N, IOK) 

Calling  Arguments: 

Name  (Attributes) 

A (real/array/ input) 


S (real/array/output) 

N (integer/ input) 

IOK  (integer/output) 


COtMON  Areas:  None 

FORTRAN  Data  Files:  None 

Requ i red  Sub rout i nes : None 

Method:  Single  precision  Cholesky  matrix  decomposition  is  performed, 

yielding  the  lower  triangular  factor  matrix  [S]  in  the  equation: 

[A]  = [S] [ST] . 

Remarks : 

(1)  This  subroutine  can  be  used  with  subroutine  BCKSUB  to  solve 
the  equation  [A]X  = B.  DECCMP  is  called  to  factor  [A]  which  permits  the 
problem  to  be  rewritten  [S] [S^]X  = B.  BCKSUB  is  called  to  solve  [S]Y  = B 
for  Y and  then  to  solve  [S^lX  = Y for  X. 

(2)  [A]  must  be  a real,  symmetric  matrix. 

(3)  The  arrays  A and  S may  overlap  whenever  A is  not  required  for 
subsequent  calculations. 


Contents 

Lower  triangle  of  matrix  to  be 
decomposed;  stored  A(1 ,1) ,A(2,1) , 

A( 2 ,2) ,A(3 ,1) A(N,N-1) ,A(N,N) . 

Length  of  array  is  ((N+l)N)/2. 
Lower  triangle  of  Cholesky  factor 
matrix;  same  order  of  storage  and 
length  as  A. 

Order  of  matrix  A. 

Flag  word.  A value  of  -1 
indicates  that  decomposition  was 
terminated  because  the  matrix  A 
was  singular. 


! 
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(4)  The  form  of  the  factor  matrix  [S]  is  compatible  with  the 
requirements  of  subroutine  BCKSUB. 
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MKMARQ  - Subroutine  Description 
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find  the  minimum  of  the  sum  of  squares  of  M functions 
MKMARQ 

CALL  MKMARQ (F , CPS ,NSIG ,M ,N ,C ,X ,FX , ITMAX , WA , I ER) 


Function:  To 
of  N variables. 

Entry  Point : 

Calling  Sequence: 

Calling  Arguments: 

Name  (Attributes) 

F (real/functior./extemal) 


EPS  (real /input) 

NSIG  (integer/input) 


M (integer/ input) 

N (integer/ input) 

C (real/input) 

X (real/array/input/output) 


Contents 

Function  to  compute  functions  to 
be  minimized,  which  will  be  called 
with  two  arguments,  F (X , I ) , where 
X is  the  vector  of  N variables  and 
I is  an  index  to  select  one  of  the 
M functions. 

First  stopping  criterion  for 
iterative  method.  If  F(X,I)<F.PS 
for  all  I,  convergence  is  assumed. 
Second  stopping  criterion  for 
iterative  method.  If  for  two 
successive  iterations  the  vector 
of  variables  X is  the  same  to  NSIG 
significant  digits,  convergence  is 
assumed. 

Number  of  functions  to  be  evaluated. 
Number  of  independent  variables  in 
vector  X. 

Algorithm  modifying  parameter; 
usually  set  to  1.0.  See  Brown ^ 
for  details. 

Vector  of  variables;  length  N. 
Initial  values  for  starting 
iteration  must  be  stored  in  this 
array  by  calling  routine.  Contains 
values  of  variables  which  mini- 
mize functions  on  output. 
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FX  (real/array/output') 


Vector  containing  the  values  of 
the  M functions  evaluated  for  the 
final  set  of  variables  in  X. 

ITMAX  (integer/ input/output)  Third  stopping  criterion  and 

iteration  counter.  If  more  than 
ITMAX  iterations  are  required  for 
convergence,  the  subroutine 
assumes  that  convergence  can  not 
be  achieved.  On  output,  ITMAX 
contains  the  number  of  iterations 
actually  performed. 

WA  (array/scratch)  A working  storage  area  KORF.  words 

long,  where 

KORF  = N (N+l) +6N+N*M+3M. 

IER  (integer/output)  Flag  word.  A value  of  zero 

indicates  that  subroutine  has 
stopped  by  meeting  first  stopping 
criterion.  A value  of  one 
indicates  that  second  stopping 
criterion  has  been  met.  If 
greater  than  one,  subroutine  has 
fai led. 

OCX  MON  Areas:  None 

FORTRAN  Data  Files:  None 

Required  Subroutines:  DECCMP , BCKSUB 

— 

Method : This  is  an  implementation  of  Brown's  Algorithm'  which  has 

been  tailored  for  use  with  the  B-spline  intersection  subroutines. 

Remarks : 

(1)  Poor  results  have  been  observed  when  this  routine  has  been 
initiated  with  zero  or  near  zero  independent  variable  vector. 

(2)  This  implementation  has  evolved  from  a similar  program  called 
ZXMARQ  which  was  marketed  at  one  time  by  International  Mathematical  and 
Statistical  Libraries,  Inc.  of  Houston,  Texas  (IMSL).  In  its  present  form 
MKMARQ  differs  from  the  IMSL  program  in  a number  of  ways  and  although  the 
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To  compute  the  matrix  product  C = A x B for  matrices 
multidimensional  elements. 


CALL  MULMAT (A , B , C , I AT , I BT , KAD , KAR , KAC , I AD , KBD , 
KBR , KBC , I BD , KCD , KCR , KCC , ICD) 


MULMAT  - Subroutine  Description 

Function: 
which  may  have 

Entry  Point : MULMAT 

Calling  Sequence: 

Calling  Arguments: 

Name  (Attributes) 

A (real/array/ input) 

B (real/array/input) 

C (real/array/output) 

IAT  (integer/ input) 

IBT  (integer/input) 

KAD, KAR, KAC  (integer/ 
input) 

IAD  (integer/input) 

KBD,KBR,KBC  (integer/ 
input) 

IBD  (integer/ input) 

KCD, KCR, KCC  (integer/ 
input i 

ICD  ( integer/ input ) 


Contents 

First  matrix  factor.  Dimensioned 
A(KAD,KAR,KAC) . 

Second  matrix  factor.  Dimensioned 
B(KBD,KBR,KBC) . 

Product  matrix.  Dimensioned  C(KCD,KCR, 
KCC). 

Flag  word.  A non-zero  value  indicates 
the  matrix  A is  to  be  transposed  prior 
to  multiplication. 

Flag  word.  A non-zero  value  indicates 
that  the  matrix  B is  to  be  transposed 
prior  to  multiplication. 

Dimensions  of  the  array  A. 

Index  to  select  the  component  of  the 
multidimensional  elements  of  array  A. 
Dimensions  of  the  array  B. 

Index  to  select  the  component  of  the 
multidimensional  elements  of  the 
array  B. 

Dimensions  of  the  array  C. 

Index  indicating  the  component  of  the 
multidimensional  elements  of  the  array 
C. 
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CC^MON  Areas : None 

FORTRAN  Data  Files:  *6  (output/ formatted) 

Required  Subroutines:  None 

Method:  Use  formula,  c. . = Z a.,b,  selecting  components  and 

1J  all  k 1K  KJ 

transposing  as  indicated. 

Remarks : 

(1)  Subroutine  will  stop  program  with  message  if  matrices  are 
not  conformable. 


N0BSF3  - Subroutine  Description 

Function:  Dummy  subroutine  to  satisfy  unused  external  references  in 

subroutine  INT3S  and  to  prevent  loading  of  subroutine  BSEVL3  and  function 
subroutine  RF3S01. 

F.ntry  Point : N0BSF1 

Calling  Sequence:  CALL  N0BSF1 


Ijytry  Point:  RF3S01 

Calling  Sequence:  CALL  RF3S01 

Calling  Arguments:  None 

Fntry  Point : BSFVL3 

Calling  Sequence:  CALL  BSFVL3 

Call ing  Arguments : None 

Entrv  Point:  RF3S03 

j. 

Calling  Sequence:  CALL  RF3S03 


Entry  Point:  CALL  RF3S04 

Calling  Sequence:  CALL  RF3S04 


ts : None 


COMMON  Areas:  None 

FORTRAN  Data  Files:  “6  (output/formatted) 

Subroutines  Required:  None 

Method:  When  this  subroutine  is  explicitly  loaded  with  subroutine 

INT3SX , the  library  loading  of  the  unused  subroutines  BSEVL3  and  RF3S01 
will  not  be  performed --conserving  memory. 

Remarks : 

rl)  Subroutine  will  stop  program  with  message  if  entry  points 
BSFVL3  and  RF3S01  are  actually  called. 
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SAMPLE  APPLICATION 


PLOTTING  B-SPLINE  SURFACES 


Although  we,  as  the  developers,  feel  that  subroutines  included  in  the 
G-PRIMF.  Basic  B-Spline  Library  are  easy  to  use,  an  example  of  a typical 
application  of  this  capability  may  be  helpful  in  establishing  a frame  of 
reference  for  the  potential  user.  The  subroutine  PLT3DS,  listed  in 
Figure  5,  draws  a plot  of  a B-spline  surface.  This  subroutine  is  from 
the  program  BHULL,  written  for  ship  hull  design,  which  is  described  in  a 
forthcoming  report.  (BHULL  also  contains  examples  of  data  fitting  and 
simple  surface  and  volume  integration  using  B-spline  functions.) 

PLT3DS  graphically  represents  a B-spline  surface  giving  the  user  the 
option  of  (1)  having  the  boundary  curves  displayed  (values  of  zero  and  one 
for  the  parameters  s and  t) , (2)  having  curves  of  constant  values  of 
either  parameter  displayed,  or  (3)  having  curves  of  constant  values  of  one 
parameter  displayed  and  then  having  curves  of  constant  values  of  the  other 
parameter  displayed  as  well. 

The  calling  sequence  for  PLT3DS  is: 

CALL  PLT3DS (A, B ,C ,MD ,M1 ,N1 , ISPERD , ITPERD , IFS , JFT , I CLEAR) 
where  the  calling  arguments  are  defined  as  follows: 


Name  (Attributes) 

A (real/array/input) 

B (array/scratch) 

C (array/scratch) 

MD,M1,N1  (integer/ input) 
ITYPE  (integer/ input) 

IFS  (integer/ input) 


Contents 

Mesh  defining  a B-spline  surface. 

Dimens ioned  A (MD ,M1 , NT) . 

Working  storage  for  PLT3DS.  Dimensioned 
B(MD,KORE)  where  KORE=MAXO (Ml ,N1) . 
Working  storage  for  subroutine  CVPLT. 
Dimensioned  C(2,K0RE+1). 

Dimensions  of  array  A. 

Code  word  indicating  type  of  B-spline 
surface  (see  Table  1) . 

Flag  word.  If  set  to  one,  curves  are 
drawn  by  varying  the  parameter  s for  a 
uniformly  spaced  sequence  of  values  of 
the  parameter  t.  If  set  to  zero. 
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SUBROUTINE  Pi  T30SI  A,B,C,NO,N  1 .HI,  XTYPC.IfS,  Jf  T, I CLEAR  I 

c 

C PI  OT  A 3-0  surface  ... 

C SURFACE  OUTLINES  ARE  AlNAYS  PLOTTEO.  IFS  AHO  Jf  T CONTROL 

C THE  PLOTTING  OF  T ANO  S CURVCS  HE  SPEC  T I VEl T . 

C 

OINENSION  A(H0,H1,N1I,  8(HD,i>,  C I 2,  1 ) 

C 

NUIP  * 6 
ILIN  * 0. 

<SCl  * 1. 

c 

C PLOT  BOUNDARIES... 

C 

IF  (PI  .IT.  1 .OR.  HI  .IT.  II  GO  TO  100 

HPT  * HULP  • Ni 

OELN  « 1.  / FlO»  TIHPT) 

HP  T * HP  T *1 
HPT  ■ HULP  • Nl 
OELN  * 1.  / FLOA  T(HPTI 
HPT  » HPT  ♦ l 
JCLEAR  * ICLEAR 


IF 

(HI  . 

EQ. 

11 

HPT 

■ 

1 

IF 

(HI  . 

CQ. 

1 1 

NPT 

a 

1 

IF 

(HI  . 

EQ. 

11 

GO 

TO 

SO 

IF 

(HI  . 

EQ. 

1) 

GO 

TO 

22 

IF 

(IFS 

.EQ. 

. 11 

GO 

TO 

11 

00  10  Ml ,2 

s * -oa« 

DO  5 1*1,  HPT 
S * S ♦ OELN 

5 CALL  BSEVll  (A.HO.Hl  ,H1,S  ,Fi.  OATU-1)  , ITTPE.Bi  1,  I I > 
CALL  C FPL  T <B, HO, HPT, Cl 
1C  JCLEAR  * 0 

11  IF  (JFT  .EQ.  II  GO  TO  22 
00  2 0 Ml, 2 
T * - OELN 
00  IS  1*1  ,HPT 
T * T ♦ OELN 

IS  CALL  BSE  VI 1 (A, HO, HI, HI, FLOAT  (K-l  I ,T  , I TY  PE  , B ( 1 , 1 1 I 
CALL  C FPL  T IB,H0,NPT,C1 
20  JCLEAR  * 0 

22  IF  (IFS  • EQ.  Cl  GO  TO  SO 

C 

c 

C PLOT  T-CURVES  ... 

C 

ISICP  * 0 
T * -OELN 
00  3 0 Mt,NPT 
T * T ♦ OELH 

IF  (ISICP  .Ne.  01  GO  To  30 
S * -OELH 
00  25  1*1 »HPT 
S * S ♦ OELH 

25  CALL  BSC  VI  1 ( A , *0,  HI  , HI , S , T,  I T TPE  , B(  1 , 1 1 I 
CALL  CVPlT  (B. NO, HPT, Cl 
XLEAP  * 0 

30  IS  HP  * HODCISKP* 1, 31 
SC  IF  (JFT  .EQ.  01  GO  TO  100 
IF  (HI  .EQ.  II  GO  TO  100 
C 
C 

C PLOT  S -CURVE S • • • 

C 

IS  KP  * 0 

s * -oan 

00  7 0 Ml, HPT 
S * S ♦ OELH 

IF  ( I SKP  .NE.  01  GO  TO  7 0 

t • -oan 

oo  so  1*1, HPT 
T * T ♦ Oan 

6 C CALL  BSE VL l (A,H0,H1,H1,S,T,  ITTPE ,81 1 ,11  1 
CALL  CVPIT  IB, HO, HPT, C) 

."•LEAR  * C 

7 C IS  HP  = HOOUSKP*  1,  11 
IOC  COHTINUE 
Rf  TURN 

C 

ENO 


l-igure  5 - Subroutine  PLT3DS 


Name  (Attributes) 


Contents 


only  the  curves  t=0  and  t=l  are  drawn. 

JFT  ( integer/ input ) Flag  word.  If  set  to  one,  curves  are 

drawn  by  varying  the  parameter  t for 
a uniformly  spaced  sequence  of  values 
for  the  parameter  s.  If  set  to  zero, 
only  the  curves  s=0  and  s=l  are  drawn. 

ICLEAR  (not  used) 

Curves  are  drawn  as  polygons  with  6*M1+1  or  6*N1+1  points  which  lie 
on  the  given  B-spline  surface.  The  four  major  processing  steps  of  this 
subroutine  are: 

(1)  The  'TO  10  loop”  draws  curves  for  t=0  and  t=l . 

(2)  The  "DO  20  loop”  draws  curves  for  s=0  and  s=l. 

(3)  The  "DO  30  loop”  draws  curves  for  t=l ./(6*N1  + 1) , 

3./ (6*N1+1) , . . . ,1 . 

(4)  The  "TO  70  loop”  draws  curves  for  s=l ./(6*M1+1) , 

3 . / (6*M1+1) , . . . ,1. 

Within  each  of  the  above  steps,  the  evaluation  subroutine  BSFVL1  is 
called  repeatedly  to  build  the  "curve  polygon"  that  is  to  be  plotted.  That 
polygon,  stored  in  the  array  B,  is  then  passed  to  subroutine  CVTLT  for 
display. 

Figure  6 shows  three  views  of  a B-spline  surface  as  drawn  by  PI.T3DS 
with  IFS=1  and  .JFT=0;  Figure  7 shows  the  same  surface  drawn  with  IFS=1 
and  JFT=1;  and  in  Figures  8 and  9 the  outlines  of  the  surfaces  have  been 
drawn  by  PLT3DS  with  IFS=0  and  JFT=0.  The  intersection  curves  and  various 
symbols  have  been  added  by  other  subroutines. 
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Figure  8 - B-Spline  Surface  Outlines  Drawn  by 
Subroutine  PLT3DS  with  IFS=0  and  JFT=0 
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Drawn  by  Subroutine  PLT3DS 


RECOMMENDED  EXTENSIONS 


This  report  represents  one  step  toward  the  creation  of  a library  of 
B-spline  mathematical  subroutines  that  will  satisfy  our  requirements  and 
the  needs  of  the  naval  scientific  community.  Our  work  in  this  field  is 
continually  uncovering  improved  methods  and  new  approaches  to  the 
problems  we  have  been  dealing  with  and  there  is  always  a list  of  modifi- 
cations and  additions  to  be  made.  Also,  for  lack  of  time,  we  have  passed 
over  some  of  the  more  general  capabilities  for  which  we  have  no  current 
applications.  This  section  lists  possible  additions  to  a future  edition 
of  the  library. 

Additional  evaluation  capabilities  should  include: 

• Calculation  of  derivatives  of  the  B-spline  functions.  This 
would  permit  the  use  of  more  efficient  numerical  integration  and  nonlinear 
minimization  algorithms. 

• Addition  of  a BSEVL2-like  subroutine  which  will  return  a set  of 
points  along  a curve  rather  than  just  one  point  per  call.  This  type  of 
subroutine  could  reduce  the  computation  time  required  to  display  curves 
and  surfaces. 

• /In  alternate  surface  definition  using  a mesh  of  triangles  as 
well  as  the  current  rectangular  mesh  definition.  Local  surface  changes 
are  more  easily  accomplished  with  a surface  defined  by  a triangular  mesh. 

Additional  fitting  capabilities  should  include: 

• A fitting  capability  which  admits  constraints  on  selected 
points  or  derivatives.  This  feature  appears  to  be  necesstn  f r fitting 
data  with  cusps  and  for  fitting  data  that  must  mesh  exactly  along  certain 
boundaries. 

• .An  option  to  use  a uniform  B-spline  basis  for  the  fitting  of 
open  curves  and  surfaces.  The  current  practice  for  open  curves  and 
surfaces  involves  the  use  of  a uniform  basis  which  is  modified  to  force 
the  B-spline  function  to  interpolate  the  ends  of  the  defining  polygon  and 
the  comers  of  a surface  mesh.  This  practice  also  constrains  the  B-spline 
function  to  have  zero  curvature  at  those  extremes.  The  uniform  basis, 
without  modification,  will  permit  arbitrary  curvature  at  the  ends,  which 
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may  be  more  appropriate  for  fitting  certain  data. 

• An  option  to  save  the  least-squares  transformation  matrices 
from  one  fitting  application  so  that  they  would  not  have  to  be  recalculated 
for  a subsequent  application.  Although  it  is  evident  from  Equations  (7) 
and  (11)  that  the  matrices  [U]  and  [W]  need  to  be  computed  only  once  for 
all  problems  of  the  same  size,  the  current  implementation  does  not  take 
advantage  of  this  property. 

• An  option  to  have  the  fitting  subroutines  compute  the  RMS 
error  in  the  fit  of  a B-spline  function  to  a given  set  of  data. 

Additional  intersection  capabilities  should  include: 

• Calculation  of  the  points  of  intersection  of  two  curves,  and 
the  intersection  of  a curve  and  a surface.  These  calculations  are 
specializations  of  the  procedures  used  for  finding  the  point  of  inter- 
section of  three  surfaces. 

• Enhancement  of  the  procedures  for  finding  the  curve  of  inter- 
section of  two  surfaces.  The  current  two-surface  intersection  capability 
is  not  as  comprehensive  as  we  require.  We  hope  to  remove  most  of  the 
restrictions  currently  imposed  and  to  allow  the  user  to  select  among 
multiple  curves  of  intersection  when  they  occur. 

• A separate  locating  capability.  The  intersection  subroutines 
now  contain  code  which  comj  s the  location  of  the  point  on  a B-spline 
curve  or  surface  that  is  closest  to  a given  spatial  point.  This  locating 
capability  should  be  made  directly  accessible  to  the  user. 

• Modification  of  the  nonlinear  solver  to  use  derivatives  of  the 
B-spline  functions  which  have  been  calculated  by  the  evaluation  subroutines 
This  should  improve  the  overall  efficiency  of  the  nonlinear  solution 
process,  since  greater  accuracy  can  be  obtained  with  no  increase  in  the 
amount  of  computation. 

• Include  curve  and  surface  operands  as  calling  arguments  of  the 
intersection  subroutines  rather  than  using  CCM10N  block  /BOPRXX/  for  the 
passing  of  operands.  This  practice,  which  is  a hold-over  from  earlier 
iterative  solution  procedures,  is  sometimes  confusing  to  the  programmer 
and  can  now  be  eliminated. 
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